!---------------------------------------------------------------------------------
! this file came from wa17 and modified by Fabrizio 07-02-2004
! standard gw_drag with modification (6) of latitude profile of gw spectrum
!---------------------------------------------------------------------------------

module gw_drag

!---------------------------------------------------------------------------------
! Purpose:
!
! Module to compute the forcing due to parameterized gravity waves. Both an 
! orographic and an internal source spectrum are considered.
!
! Author: Byron Boville
!
!---------------------------------------------------------------------------------
  use shr_kind_mod,   only : r8 => shr_kind_r8
  use ppgrid,         only : pcols, pver, pverp
  use constituents,   only : pcnst
  use physics_types,  only : physics_state, physics_ptend
  use spmd_utils,     only : masterproc
  use cam_history,    only : outfld
  use physconst,      only : gravit, pi
  use cam_logfile,    only : iulog
  use abortutils,     only : endrun
 
  implicit none

  save
  private                         ! Make default type private to the module
!
! PUBLIC: interfaces
!
  public :: gw_inti                  ! Initialization
  public :: gw_intr                  ! interface to actual parameterization
  public :: gw_drag_readnl           ! Read namelist
  public :: gw_drag_register
  public :: idx_zmdt

!
! PRIVATE: Rest of the data and interfaces are private to this module
!
  integer,  parameter :: pgwv = 32           ! number of waves allowed (one-signed only)
  real(r8), parameter :: dc   = 2.5_r8       ! bin width for spectrum
  real(r8), parameter :: dc2  = 5.0_r8       ! bin width for output spectrum (for outputting tauxs)
  real(r8), parameter :: dca  = 0.1_r8       ! integration interval to get bin average
  real(r8), parameter :: spec_fac=10._r8/dc  ! number of bins compared to dc=10	!bab
  real(r8), parameter :: c0   = 30._r8       ! width of gaussian in phase speed
  real(r8) :: fav(-pgwv:pgwv)                ! average value of gaussian over bin
  integer :: kbotbg, kbotoro                 ! interface of gwd source
  integer :: ktopbg, ktoporo                 ! top interface of gwd region


!+Mods for Beres04
  integer :: bkbotbg                         ! interface of gwd beres source
  integer :: k700                            ! 700 mb index (used in the gw_bgnd_beres)
  real(r8), parameter :: effgw_beres=0.1_r8  ! beres source tendency efficiency
!-Mods for Beres04

!+Mods for CM
  real(r8) :: cref(-pgwv:pgwv)               ! Reference phase speed spectrum
  real(r8) :: cref2(-pgwv:pgwv)               ! Reference phase speed spectrum 2 (just for outputting tauxs)
  real(r8) :: frontgfc                       ! frontogenesis function critical threshold
!-Mods for CM


  real(r8) :: alpha(0:pver)                  ! newtonian cooling coefficients
  real(r8) :: cpair                          ! specific heat of dry air (constant p)
  real(r8) :: annphase                       ! phase of annual cycle for taubgnd (day of max)
  real(r8) :: annavg_n, annamp_n             ! NH annual average and amplitude for taubgnd
  real(r8) :: annavg_s, annamp_s             ! SH annual average and amplitude for taubgnd
  real(r8) :: dback                          ! background diffusivity
  real(r8) :: effkwv                         ! effective wavenumber (fcrit2*kwv)
  real(r8) :: effgw_spec                     ! efficiency factor for spectrum !bab
  real(r8) :: effgw_oro                      ! efficiency factor for orographic wave
  real(r8), parameter :: fracldv    = 0._r8  ! fraction of stress deposited in low level region
  real(r8) :: g                   ! acceleration of gravity
  real(r8) :: kwv                 ! effective horizontal wave number
  real(r8) :: lzmax               ! maximum vertical wavelength at source

  real(r8) :: mxasym              ! max asymmetry between tau(c) and tau(-c)
  real(r8) :: mxrange             ! max range of tau for all c

  real(r8) :: n2min               ! min value of bouyancy frequency
  real(r8) :: fcrit2              ! critical froude number
  real(r8) :: oroko2              ! 1/2 * horizontal wavenumber
  real(r8) :: orohmin             ! min surface displacment height for orographic waves
  real(r8) :: orovmin             ! min wind speed for orographic waves
  real(r8) :: r                   ! gas constant for dry air
  real(r8) :: rog                 ! r / g
  real(r8) :: taubgnd             ! background source strength (/tauscal)
  real(r8) :: taumin              ! minimum (nonzero) stress
  real(r8) :: tauscal             ! scale factor for background stress source
  real(r8) :: tndmax              ! maximum wind tendency
  real(r8) :: umcfac              ! factor to limit tendency to prevent reversing u-c
  real(r8) :: ubmc2mn             ! min (u-c)**2
  real(r8) :: zldvcon             ! constant for determining zldv from tau0

  integer, parameter :: nalph=66  ! no. of levels of pre-calculated newtonian cooling
  real(r8) :: alpha0(nalph)       ! newtonian cooling (1/day)
  real(r8) :: palph(nalph)        ! pressure levs of pre-calculated newtonian cooling (hpa)

  real(r8), parameter :: d2r  = pi/180._r8

  data alpha0 /  1.896007_r8   , 1.196965_r8   , 0.7251356_r8  ,  0.6397463_r8   , &
                 0.5777858_r8  , 0.5712274_r8  , 0.6836302_r8  ,  0.6678557_r8   , &
                 0.5683219_r8  , 0.4754283_r8  , 0.3960519_r8  ,  0.332022_r8    , &
                 0.2497581_r8  , 0.168667_r8   , 0.1323903_r8  ,  0.1257139_r8   , &
                 0.1069889_r8  , 0.09873954_r8 , 0.09215571_r8 ,  0.09398635_r8  , &
                 0.1061087_r8  , 0.1294598_r8  , 0.1544743_r8  ,  0.1648226_r8   , &
                 0.1687332_r8  , 0.1691513_r8  , 0.1664987_r8  ,  0.159048_r8    , &
                 0.149292_r8   , 0.1351563_r8  , 0.1174998_r8  ,  0.09913579_r8  , &
                 0.08300615_r8 , 0.0707_r8     , 0.0615588_r8  ,  0.0542623_r8   , &
                 0.0478562_r8  , 0.04132157_r8 , 0.03454087_r8 ,  0.02296682_r8  , &
                 0.006723819_r8, 0.02164464_r8 , 0.05756261_r8 ,  0.003844868_r8 , &
                 0.02929285_r8 , 0.006627098_r8, 0.04558291_r8 ,  0.02042176_r8  , &
                 0.00000000_r8 , 0.005880283_r8, 0.00689498_r8 ,  0.01343466_r8  , &
                 0.00000000_r8 , 0.03415992_r8 , 0.02855049_r8 ,  0.01688839_r8  , &
                 0.0272628_r8  , 0.02772121_r8 , 0.02135626_r8 ,  0.04863235_r8  , &
                 0.04568304_r8 , 0.00000000_r8 , 0.009604108_r8,  0.00000000_r8  , &
                 0.00000000_r8 , 0.00000000_r8 /

  data palph /   5.11075e-06_r8 , 9.8269e-06_r8   , 1.620185e-05_r8 , 2.671225e-05_r8 , & 
                 4.4041e-05_r8  , 7.261275e-05_r8 , 0.000119719_r8  , 0.00019738_r8   , &
                 0.0003254225_r8, 0.0005365325_r8 , 0.0008846025_r8 , 0.001458458_r8  , &
                 0.002404575_r8 , 0.00397825_r8   , 0.006556825_r8  , 0.01081382_r8   , &
                 0.017898_r8    , 0.02955775_r8   , 0.04873075_r8   , 0.07991075_r8   , &
                 0.1282732_r8   , 0.19812_r8      , 0.292025_r8     , 0.4101675_r8    , & 
                 0.55347_r8     , 0.73048_r8      , 0.9559475_r8    , 1.244795_r8     , &  
                 1.61285_r8     , 2.079325_r8     , 2.667425_r8     , 3.404875_r8     , &
                 4.324575_r8    , 5.4654_r8       , 6.87285_r8      , 8.599725_r8     , & 
                 10.70705_r8    , 13.26475_r8     , 16.35175_r8     , 20.05675_r8     , &
                 24.479_r8      , 29.728_r8       , 35.92325_r8     , 43.19375_r8     , &
                 51.6775_r8     , 61.5205_r8      , 72.8745_r8      , 85.65715_r8     , &
                 100.5147_r8    , 118.2503_r8     , 139.1154_r8     , 163.6621_r8     , &
                 192.5399_r8    , 226.5132_r8     , 266.4812_r8     , 313.5013_r8     , &
                 368.818_r8     , 433.8952_r8     , 510.4553_r8     , 600.5242_r8     , &
                 696.7963_r8    , 787.7021_r8     , 867.1607_r8     , 929.6489_r8     , & 
                 970.5548_r8    , 992.5561_r8 /

!+ Mods for Beres04
! BERES parameterization arrays
  integer, parameter :: MAXH = 20                 ! max values of heating depth
  integer, parameter :: MAXUH = 40                ! max value of mean wind in heating
  real(r8) :: mfcc(MAXH,-MAXUH:MAXUH,-PGWV:PGWV)  ! stores source spectra
  character(len=256) :: gw_drag_file = 'Beres04_file'
!-Mods for Beres04

  integer :: idx_zmdt

contains

!===============================================================================

  subroutine gw_drag_register
    use ppgrid,       only : pver
    use phys_buffer,  only : pbuf_add

    ! add fields to physics buffer
    call pbuf_add('ZMDT', 'physpkg',1,pver,1,idx_zmdt)
  end subroutine gw_drag_register
  
  subroutine gw_drag_readnl(nlfile)

    use namelist_utils,  only: find_group_name
    use units,           only: getunit, freeunit
    use mpishorthand

    character(len=*), intent(in) :: nlfile  ! filepath for file containing namelist input

    ! Local variables
    integer :: unitn, ierr
    character(len=*), parameter :: subname = 'gw_drag_readnl'

    namelist /gw_drag_nl/ gw_drag_file
    !-----------------------------------------------------------------------------

    if (masterproc) then
       unitn = getunit()
       open( unitn, file=trim(nlfile), status='old' )
       call find_group_name(unitn, 'gw_drag_nl', status=ierr)
       if (ierr == 0) then
          read(unitn, gw_drag_nl, iostat=ierr)
          if (ierr /= 0) then
             call endrun(subname // ':: ERROR reading namelist')
          end if
       end if
       close(unitn)
       call freeunit(unitn)
    end if

#ifdef SPMD
    ! Broadcast namelist variables
    call mpibcast (gw_drag_file, len(gw_drag_file), mpichar, 0, mpicom)
#endif

  end subroutine gw_drag_readnl

  subroutine gw_inti (cpairx, cpwv, gx, rx, hypi)
!-----------------------------------------------------------------------
! Time independent initialization for multiple gravity wave parameterization.
!-----------------------------------------------------------------------
    use cam_history,      only : addfld, add_default, phys_decomp
    use dycore,           only : get_resolution
    use interpolate_data, only : lininterp
    use phys_buffer,      only : pbuf
    use infnan

!------------------------------Arguments--------------------------------
    real(r8), intent(in) :: cpairx                ! specific heat of dry air (constant p)
    real(r8), intent(in) :: cpwv                  ! specific heat of water vapor (constant p)
    real(r8), intent(in) :: gx                    ! acceleration of gravity
    real(r8), intent(in) :: rx                    ! gas constant for dry air
    real(r8), intent(in) :: hypi(pver+1)          ! reference interface pressures

!---------------------------Local storage-------------------------------
    integer :: k
    integer :: m, n, i, l
    character(len=80) fnum, fname, dumc1, dumc2, dumc1x, dumc1y   ! dummies
    real(r8) :: ccn
    real(r8) :: cmn,cmx                           ! bin edges (max and min values of c)

!-----------------------------------------------------------------------
! Copy model constants
    cpair  = cpairx
    g      = gx
    r      = rx

! Set MGWD constants
    kwv    = 6.28e-5_r8       ! 100 km wave length
    dback  = 0.05_r8          ! background diffusivity
    fcrit2 = 1.0_r8           ! critical froude number squared		
    tauscal= 0.001_r8         ! scale factor for background stress

    select case (get_resolution())
    case ('4x5')
       effgw_oro  = 0.125_r8
       effgw_spec = 1.0_r8	! rrg
       taubgnd    = 1.0_r8	! jhr
       frontgfc   = 7.5e-16_r8
    case ('1.9x2.5')
       effgw_oro  = 0.125_r8
!      effgw_spec = 1.0_r8 / spec_fac 
       effgw_spec = 1.0_r8	! fs (04/27/2008)
       taubgnd    = 1.5_r8
       frontgfc   = 1.25E-15_r8  !jhr
    case default
       effgw_oro  = 0.125_r8
       effgw_spec = 1.0_r8
       taubgnd    = 1.5_r8
       frontgfc   = 1.25E-15_r8  !jhr
    end select

    if (masterproc) then
       write(iulog,*) ' '
       write(iulog,*) 'gw_inti: gw spectrum taubgnd, effgw_spec = ', taubgnd, effgw_spec
       write(iulog,*) ' '
    end if

    zldvcon = 10._r8           ! constant for determining zldv
    lzmax   = 7.e3_r8          ! maximum vertical wavelength at source (m)

!+Mods for Beres04
    pbuf(idx_zmdt) %fld_ptr = nan
!-Mods for Beres04
!+Mods for CM
    cref(:) = 0._r8
    cref2(:) = 0._r8
!-Mods for CM

! the following block by BAB
!+++ bab

! If a spectrum is being computed, integrate gaussian over bins

    if (pgwv > 0) then
       write(iulog,*) 'GW_DRAG: pgwv = ', pgwv				! rrg
       do l = -pgwv, pgwv
          cref(l) = dc * l                                              ! fs
          cref2(l) = dc2 * l                                            ! jhr	
! Lower bound of bin
          cmn = cref(l) - 0.5_r8*dc
          cmx = cref(l) + 0.5_r8*dc
! Loop over integration intervals in bin
          fav(l) = 0.5_r8 * dca * (exp(-(cmn/c0)**2) + exp(-(cmx/c0)**2))
          do n = 1,nint(dc/dca)-1
             fav(l) = fav(l) + dca * exp(-((cmn+n*dca)/c0)**2)
          enddo
          fav(l) = fav(l) / dc
          if (masterproc) then
             write (iulog,*) 'GW_DRAG: c(l), fav(l)= ', cref(l), fav(l)	! rrg
          end if
       enddo
    else
       fav(0) = 1._r8
    end if

!--- bab

! pre-calculated newtonian damping: 
!     * convert to 1/s
!     * ensure it is not smaller than 1e-6
!     * convert palph from hpa to pa

  do k=1,nalph
     alpha0(k) = alpha0(k) / 86400._r8
     if (alpha0(k) .lt. 1.e-6_r8) alpha0(k) = 1.e-6_r8
     palph(k) = palph(k) *1.e2_r8
  enddo

! interpolate to current vertical grid and obtain alpha

  call lininterp (alpha0  ,palph, nalph , alpha  , hypi , pverp)
  if (masterproc) then
     write (iulog,*) 'gw_inti: newtonian damping (1/day):'
     write (iulog,fmt='(a4,a12,a10)') ' k  ','  hypi      ','  alpha   '
     do k=0,pver
        write (iulog,fmt='(i4,1e12.5,1f10.2)') k,hypi(k+1),alpha(k)*86400.
     enddo
  endif

! set radiative damping times
!!$    do k = 0, pver
!!$       alpha(k) = 1.e-6       ! about 10 days.
!!$    end do

! Min and max values to keep things reasonable
    mxasym = 0._r8		     ! max factor of 10 from |tau(c)| to |tau(-c)|
    mxrange= 0._r8                   ! factor of 100 from max to min |tau(c)|
    n2min   = 1.e-8_r8               ! min value of Brunt-Vaisalla freq squared
    orohmin = 10._r8                 ! min surface displacement for orographic wave drag
    orovmin = 2._r8                  ! min wind speed for orographic wave drag
    taumin  = 1.e-10_r8              ! min stress considered > 0
    tndmax  = 400._r8 / 86400._r8    ! rrg: max permitted tendency (500 m/s/day)
    umcfac  = 0.5_r8                 ! max permitted reduction in u-c
    ubmc2mn = 0.01_r8                ! min value of (u-c)^2

! Determine other derived constants
    oroko2 = 0.5_r8 * kwv
    effkwv = fcrit2 * kwv
    rog    = r/g

! Determine the bounds of the background and orographic stress regions
    ktopbg  = 0
    kbotoro = pver
    do k = 0, pver
       if (hypi(k+1) .lt. 50000._r8) kbotbg  = k    ! spectrum source at 500 mb
       if (hypi(k+1) .lt. 50000._r8) kbotbg  = k    ! spectrum source at 500 mb
!+ Mods for Beres04
       if (hypi(k+1) .lt. 10000.) bkbotbg  = k+1    ! beres spectrum source at 100 mb
       if (hypi(k+1) .lt. 70000.) k700    = k+1     ! 700 mb index
!- Mods for Beres04
    end do
    ktoporo = 0

    if (masterproc) then
       write (iulog,*) 'KTOPBG  =',ktopbg
       write (iulog,*) 'KBOTBG  =',kbotbg
       write (iulog,*) 'KTOPORO =',ktoporo
       write (iulog,*) 'KBOTORO =',kbotoro
!+Mods for Beres04
       write (iulog,*) 'BKBOTBG =',bkbotbg
       write (iulog,*) 'K700    =',k700
!+Mods for Beres04
    end if

! Declare history variables for orgraphic term
    call addfld ('TTGW','K/s     ',pver, 'A','T tendency - gravity wave drag',phys_decomp)
    call addfld ('UTGWORO ','m/s2    ',pver, 'A','U tendency - orographic gravity wave drag',phys_decomp)
    call addfld ('VTGWORO ','m/s2    ',pver, 'A','V tendency - orographic gravity wave drag',phys_decomp)
    call addfld ('TAUGWX  ','N/m2    ',1,    'A','Zonal gravity wave surface stress',        phys_decomp)
    call addfld ('TAUGWY  ','N/m2    ',1,    'A','Meridional gravity wave surface stress',   phys_decomp)

!   call add_default ('TTGW', 1, ' ')
    call add_default ('UTGWORO ', 1, ' ')
    call add_default ('VTGWORO ', 1, ' ')
    call add_default ('TAUGWX  ', 1, ' ')
    call add_default ('TAUGWY  ', 1, ' ')

! Declare history variables for spectrum
    if (pgwv > 0) then
       call addfld ('UTGWSPEC','m/s2    ',pver, 'A','U tendency - gravity wave spectrum',       phys_decomp)
       call addfld ('VTGWSPEC','m/s2    ',pver, 'A','V tendency - gravity wave spectrum',       phys_decomp)
       call add_default ('UTGWSPEC', 1, ' ')
       call add_default ('VTGWSPEC', 1, ' ')

       call addfld ('BUTGWSPEC ','m/s2      ',pver, 'A','Beres U tendency',phys_decomp)
       call addfld ('BVTGWSPEC ','m/s2      ',pver, 'A','Beres V tendency',phys_decomp)
       call add_default ('BUTGWSPEC', 1, ' ')
       call add_default ('BVTGWSPEC', 1, ' ')

! Gravity Wave Tendencies Separated by wave phase speed

       call addfld ('UTEND1  ','m/s2    ',pver, 'A','U tendency   c < -40',                phys_decomp)
       call addfld ('UTEND2  ','m/s2    ',pver, 'A','U tendency  -40 < c < -15',                phys_decomp)
       call addfld ('UTEND3  ','m/s2    ',pver, 'A','U tendency  -15 < c <  15',                phys_decomp)
       call addfld ('UTEND4  ','m/s2    ',pver, 'A','U tendency   15 < c <  40',                phys_decomp)
       call addfld ('UTEND5  ','m/s2    ',pver, 'A','U tendency   40 < c ',                phys_decomp)


       call addfld ('BUTEND1  ','m/s2    ',pver, 'A','U tendency   c < -40',                phys_decomp)
       call addfld ('BUTEND2  ','m/s2    ',pver, 'A','U tendency  -40 < c < -15',                phys_decomp)
       call addfld ('BUTEND3  ','m/s2    ',pver, 'A','U tendency  -15 < c <  15',                phys_decomp)
       call addfld ('BUTEND4  ','m/s2    ',pver, 'A','U tendency   15 < c <  40',                phys_decomp)
       call addfld ('BUTEND5  ','m/s2    ',pver, 'A','U tendency   40 < c',                phys_decomp)


       call addfld ('NETDT  ','K/s   ',pver, 'A','Net heating rate',                             phys_decomp)
       call addfld ('MAXQ0  ','K/day   ',1  ,  'A','Max column heating rate',                      phys_decomp)   
       call addfld ('HDEPTH  ','km    ',1,    'A','Heating Depth',                                phys_decomp)
       call addfld ('BEMF  ','Pa    ',pver,    'A','Beres Eastward MF',                                phys_decomp)
       call addfld ('BWMF  ','Pa    ',pver,    'A','Beres Westward MF',                                phys_decomp)
       call addfld ('BNMF  ','Pa    ',pver,    'A','Beres Northward MF',                                phys_decomp)
       call addfld ('BSMF  ','Pa    ',pver,    'A','Beres Southward MF',                                phys_decomp)
       call addfld ('EMF  ','Pa    ',pver,    'A','Eastward MF',                                phys_decomp)
       call addfld ('WMF  ','Pa    ',pver,    'A','Westward MF',                                phys_decomp)
       call addfld ('NMF  ','Pa    ',pver,    'A','Northward MF',                                phys_decomp)
       call addfld ('SMF  ','Pa    ',pver,    'A','Southward MF',                                phys_decomp)

       call add_default ('UTEND1   ', 1, ' ')
       call add_default ('UTEND2   ', 1, ' ')
       call add_default ('UTEND3   ', 1, ' ')
       call add_default ('UTEND4   ', 1, ' ')
       call add_default ('UTEND5   ', 1, ' ')


       call add_default ('BUTEND1   ', 1, ' ')
       call add_default ('BUTEND2   ', 1, ' ')
       call add_default ('BUTEND3   ', 1, ' ')
       call add_default ('BUTEND4   ', 1, ' ')
       call add_default ('BUTEND5   ', 1, ' ')

       call add_default ('NETDT    ', 1, ' ')
       call add_default ('HDEPTH   ', 1, ' ')
       call add_default ('MAXQ0    ', 1, ' ')

    end if

    call addfld ('TTGWSDF' ,'K/s   ',pver, 'A','t tendency - gw spec: diffusion term',phys_decomp)
    call addfld ('TTGWSKE' ,'K/s   ',pver, 'A','t tendency - gw spec: kinetic energy conversion term',phys_decomp)
!   call add_default ('TTGWSDF', 1, ' ')
!   call add_default ('TTGWSKE', 1, ' ')

    call addfld ('EKGWSPEC' ,'M2/S   ',pverp, 'A','effective Kzz due to gw spectrum',phys_decomp)
    call add_default ('EKGWSPEC', 1, ' ')

    if (pgwv > 0) then
       call addfld ('TAUE' ,'Pa   ',pverp, 'A','CM Eastward Reynolds stress',phys_decomp)
       call addfld ('TAUW' ,'Pa   ',pverp, 'A','CM Westward Reynolds stress',phys_decomp)
       call addfld ('TAUNET' ,'Pa   ',pverp, 'A','CM E+W Reynolds stress',phys_decomp)
       call addfld ('TAUN' ,'Pa   ',pverp, 'A','CM Northward Reynolds stress',phys_decomp)
       call addfld ('TAUS' ,'Pa   ',pverp, 'A','CM Southward Reynolds stress',phys_decomp)
       call add_default ('TAUE', 1, ' ')
       call add_default ('TAUW', 1, ' ')
       call add_default ('TAUNET', 1, ' ')
       call add_default ('TAUN', 1, ' ')
       call add_default ('TAUS', 1, ' ')

       call addfld ('BTAUE' ,'Pa   ',pverp, 'A','Beres Eastward Reynolds stress',phys_decomp)
       call addfld ('BTAUW' ,'Pa   ',pverp, 'A','Beres Westward Reynolds stress',phys_decomp)
       call addfld ('BTAUN' ,'Pa   ',pverp, 'A','Beres Northward Reynolds stress',phys_decomp)
       call addfld ('BTAUS' ,'Pa   ',pverp, 'A','Beres Southward Reynolds stress',phys_decomp)
       call add_default ('BTAUE', 1, ' ')
       call add_default ('BTAUW', 1, ' ')
       call add_default ('BTAUN', 1, ' ')
       call add_default ('BTAUS', 1, ' ')

!+Mods for CM
       call addfld ('FRONTGF', 'K^2/M^2/S', pver, 'A', 'Frontogenesis function at gws src level', phys_decomp)
       call add_default ('FRONTGF', 1, ' ')
       call addfld ('FRONTGFA', 'K^2/M^2/S', pver, 'A', 'Frontogenesis function at gws src level', phys_decomp)
       call add_default ('FRONTGFA', 1, ' ')
!=Mods for CM



! Gravity wave source spectra by wave number

! Positive phase velocities
      l=-1
      do m=0,PGWV
         ccn=float(m)*dc2
         if (ccn <= cref2(pgwv)) then
            l=l+1
            write (fnum,fmt='(i3)') 100+l
            dumc1x='TAUXSp'//fnum(2:3)
            dumc1y='TAUYSp'//fnum(2:3)
            write (fnum,fmt='(f6.2,a5,f6.2,a4)') ccn,' m/s'
            dumc2='tau at c= '//fnum
	    call addfld (dumc1x(1:8),'Pa   ',pver, 'A',dumc2,phys_decomp)
            call addfld (dumc1y(1:8),'Pa   ',pver, 'A',dumc2,phys_decomp)

         end if
      end do


! Positive phase velocities
      l=-1
      do m=0,PGWV
         ccn=float(m)*dc
         if (ccn <= cref(pgwv)) then
            l=l+1
            write (fnum,fmt='(i3)') 100+l
            dumc1x='BTAUXSp'//fnum(2:3)
            dumc1y='BTAUYSp'//fnum(2:3)
            write (fnum,fmt='(f6.2,a5,f6.2,a4)') ccn,' m/s'
            dumc2='Beres tau at c= '//fnum
            call addfld (dumc1x(1:9),'Pa   ',pver, 'A',dumc2,phys_decomp)
            call addfld (dumc1y(1:9),'Pa   ',pver, 'A',dumc2,phys_decomp)

         end if
      end do


! Negative phase velocities
      l=0
      do m=-1,-PGWV,-1
         ccn=float(m)*dc2
         if (ccn >= cref2(-pgwv)) then
            l=l+1
            write (fnum,fmt='(i3)') 100+l
            dumc1x='TAUXSn'//fnum(2:3)
            dumc1y='TAUYSn'//fnum(2:3)
            write (fnum,fmt='(f6.2,a5,f6.2,a4)') ccn,' m/s'
            dumc2='tau at c= '//fnum
            call addfld (dumc1x(1:8),'Pa   ',pver, 'A',dumc2,phys_decomp)
            call addfld (dumc1y(1:8),'Pa   ',pver, 'A',dumc2,phys_decomp)

         end if
      end do



! Negative phase velocities
      l=0
      do m=-1,-PGWV,-1
         ccn=float(m)*dc
         if (ccn >= cref(-pgwv)) then
            l=l+1
            write (fnum,fmt='(i3)') 100+l
            dumc1x='BTAUXSn'//fnum(2:3)
            dumc1y='BTAUYSn'//fnum(2:3)
            write (fnum,fmt='(f6.2,a5,f6.2,a4)') ccn,' m/s'
            dumc2='Beres tau at c= '//fnum
            call addfld (dumc1x(1:9),'Pa   ',pver, 'A',dumc2,phys_decomp)
            call addfld (dumc1y(1:9),'Pa   ',pver, 'A',dumc2,phys_decomp)

         end if
      end do



    endif

!+Mods for Beres04
! Initialization of BERES' parameterization parameters
    call gw_inti_beres
!-Mods for Beres04

    return

  end  subroutine gw_inti

!+Mods for Beres04
!==============================================================================
  subroutine gw_inti_beres

   use ioFileMod,        only: getfil
#if ( defined SPMD )
   use mpishorthand
#endif

   include 'netcdf.inc'

! VARIABLES NEEDED TO READ IN TABLE OF SPECTRA FROM FILE

   integer ih, iw,ip                ! height,wind, wave look-up table indices
   integer ncid, varid, result       ! netcdf file id, variable id, result of netcdf call
   real(r8) :: tmpmfcc(20,-40:40,-40:40)  ! stores source spectra
   character(len=256) :: gw_drag_file_loc ! local filepath of gw_drag_file

   !----------------------------------------------------------------------
   ! read in look-up table for source spectra
   !-----------------------------------------------------------------------

   if (masterproc) then
      
     call getfil(gw_drag_file, gw_drag_file_loc)
     result = NF_OPEN (gw_drag_file_loc,0,ncid)

      if (result .ne. 0) then
         write(iulog,*) 'Error reading in netcdf file ',gw_drag_file,'.  ', NF_STRERROR(result)
         write(iulog,*) 'Be sure that the variable gw_drag_file is set properly in gw_drag.F90.'
         call endrun
      endif

      result = NF_INQ_VARID (ncid,'mfcc',varid)

      if (result .ne. 0) then
         write(iulog,*) 'Error reading data from ',gw_drag_file,'.  ', NF_STRERROR(result)
         call endrun
      endif

      result = NF_GET_VAR_DOUBLE (ncid,varid,tmpmfcc)
      result = NF_CLOSE (ncid)

      mfcc(:,:,:) = tmpmfcc(:,:,-PGWV:PGWV)   ! reduce the amount of wave numbers needed

      write(iulog,*) 'Read-in source spectra from file'
      write (iulog,*) 'MFCC=',maxval(mfcc),minval(mfcc)

   endif
 !     Broadcast results
#ifdef SPMD
   call mpibcast (mfcc, MAXH*(2*PGWV+1)*(2*MAXUH+1),mpir8,0,mpicom)
#endif

 end subroutine gw_inti_beres

!-Mods for Beres04

!===============================================================================

  subroutine set_phase_speeds (lchnk, ncol, ngwv, ubi , c, bflag ) 

    integer, intent(in) :: lchnk
    integer, intent(in) :: ncol
    integer, intent(in) :: ngwv

    logical, intent(in) :: bflag

    real(r8), intent(in) :: ubi(pcols,0:pver)
    
    real(r8), intent(out) :: c(pcols,-pgwv:pgwv)

!------------------------------------------------------------------------------
    integer l, i, n
    
    c(:,:) = 0._r8
    ! Check if dealing with the spectrum or the orographic component
    if ( ngwv /= 0) then

       ! Now check if dealing with the Beres04 component of the spectrum
       if (bflag) then

          do l=-ngwv,ngwv
             do i =1, ncol
                c(i,l) = cref(l)
             end do
          end do

       else

          do l=-ngwv,ngwv
             do i =1, ncol
                c(i,l) = abs(ubi(i,kbotbg)) + cref(l)
             end do
          end do

       end if

    else

       do i=1, ncol
          c(i,0) = 0._r8
       end do

    end if


  end subroutine set_phase_speeds

  

!===============================================================================

  subroutine gw_intr (state,  sgh,  pblh,   dt, ptend, landfrac, &
                      kvh )

!-----------------------------------------------------------------------
! Interface for multiple gravity wave drag parameterization.
!-----------------------------------------------------------------------
!------------------------------Arguments--------------------------------

    use time_manager, only: get_nstep,get_step_size

    real(r8), intent(in) :: sgh(pcols)            ! standard deviation of orography
    real(r8), intent(in) :: pblh(pcols)           ! planetary boundary layer height
    real(r8), intent(in) :: dt                    ! time step
    real(r8), intent(in) :: landfrac(pcols)        ! Land fraction

    type(physics_state), intent(in) :: state      ! physics state structure
    type(physics_ptend), intent(inout):: ptend    ! parameterization tendency structure

    real(r8), intent(in) :: kvh(pcols,0:pver)           ! molecular thermal diffusivity

!---------------------------Local storage-------------------------------
    integer :: lchnk                              ! chunk identifier
    integer :: ncol                               ! number of atmospheric columns

    integer :: i,k                                ! loop indexes
    integer :: kldv(pcols)                        ! top interface of low level stress divergence region
    integer :: kldvmn                             ! min value of kldv
    integer :: ksrc(pcols)                        ! index of top interface of source region
    integer :: ksrcmn                             ! min value of ksrc

    real(r8) :: ttgw(pcols,pver)                  ! temperature tendency
    real(r8) :: utgw(pcols,pver)                  ! zonal wind tendency
    real(r8) :: vtgw(pcols,pver)                  ! meridional wind tendency

    real(r8) :: ni(pcols,0:pver)                  ! interface Brunt-Vaisalla frequency
    real(r8) :: nm(pcols,pver)                    ! midpoint Brunt-Vaisalla frequency
    real(r8) :: rdpldv(pcols)                     ! 1/dp across low level divergence region
    real(r8) :: rhoi(pcols,0:pver)                ! interface density
    real(r8) :: tau(pcols,-pgwv:pgwv,0:pver)      ! wave Reynolds stress
    real(r8) :: tau0x(pcols)                      ! c=0 sfc. stress (zonal)
    real(r8) :: tau0y(pcols)                      ! c=0 sfc. stress (meridional)
    real(r8) :: ti(pcols,0:pver)                  ! interface temperature
    real(r8) :: ubi(pcols,0:pver)                 ! projection of wind at interfaces
    real(r8) :: ubm(pcols,pver)                   ! projection of wind at midpoints
    real(r8) :: xv(pcols)                         ! unit vectors of source wind (x)
    real(r8) :: yv(pcols)                         ! unit vectors of source wind (y)

    character(len=8) :: fnum,fname                ! dummy characters
    integer  ic,l,m                               ! dummy integers
    real(r8) qtgw(pcols,pver,pcnst)               ! constituents tendencies

    real(r8) :: taue(pcols,0:pver)                ! Reynolds stress for eastward propagating waves
    real(r8) :: tauw(pcols,0:pver)                ! Reynolds stress for westward propagating waves
    real(r8) :: taun(pcols,0:pver)                ! Reynolds stress for northward propagating waves
    real(r8) :: taus(pcols,0:pver)                ! Reynolds stress for southward propagating waves

    real(r8) :: btaue(pcols,0:pver)                ! Reynolds stress for eastward propagating waves Beres scheme
    real(r8) :: btauw(pcols,0:pver)                ! Reynolds stress for westward propagating waves Beres scheme
    real(r8) :: btaun(pcols,0:pver)                ! Reynolds stress for northward propagating waves Beres scheme
    real(r8) :: btaus(pcols,0:pver)                ! Reynolds stress for southward propagating waves Beres scheme

    real(r8) :: c(pcols,-pgwv:pgwv)                       ! spectrum phase speeds for each column
! C.-C. Chen
    integer :: src_level(pcols)                    ! index of gravity wave source
!-----------------------------------------------------------------------------
    lchnk = state%lchnk
    ncol  = state%ncol

    taue(:,:) = 0._r8
    tauw(:,:) = 0._r8
    taun(:,:) = 0._r8
    taus(:,:) = 0._r8

    btaue(:,:) = 0._r8
    btauw(:,:) = 0._r8
    btaun(:,:) = 0._r8
    btaus(:,:) = 0._r8

! Profiles of background state variables
    call gw_prof(lchnk, ncol, &
         state%u   , state%v   , state%t   , state%pmid   , state%pint, &
         rhoi      , ni        , ti        , nm           )



!-----------------------------------------------------------------------------
! Non-orographic backgound gravity wave spectra
!-----------------------------------------------------------------------------
    if (pgwv >0) then

!+Mods for Beres04
! Determine wave sources for Beres04 scheme
! set phase speeds for this chunk
       call set_phase_speeds (lchnk, ncol, PGWV,  ubi , c , .true.)

       call gw_bgnd_beres (lchnk          , ncol       ,                     &
            state%u    , state%v    , state%t    , state%pmid , state%pint , &
            state%pdel , state%rpdel, state%lnpint,kldv       , kldvmn     , &
            ksrc       , ksrcmn     , rdpldv     , tau        , ubi        , &
            ubm        , xv         , yv         , PGWV       , &
            state%zm, src_level)

! C.-C. Chen
!       src_level(:) = bkbotbg

! Solve for the drag profile with Beres source spectrum
       call gw_drag_prof (lchnk     , ncol       ,                     &
            PGWV       , bkbotbg     , ktopbg     , state%u    , state%v     , &
            state%t    , state%pint , state%pdel , state%rpdel, state%lnpint ,&
            rhoi       , ni         , ti         , nm         , dt          , &
            kldv       , kldvmn     , ksrc       , ksrcmn     , rdpldv      , &
            tau        , ubi        , ubm        , xv         , yv          , &
            effgw_beres , utgw       , vtgw       , tau0x      , tau0y      , &
            kvh        , state%pmid , state%q    , state%s    , ttgw       , qtgw, .true.     , &
            taue       , tauw       , taun       , taus       , c, src_level)

! Store constituents tendencies
!caf       do m=1, pcnst+pnats
       do m=1, pcnst
          do k = 1, pver
             do i = 1 , ncol
                ptend%q(i,k,m) = qtgw(i,k,m)
             end do
          end do

       end do                          ! another constituent

! add the momentum tendencies to the output tendency arrays
       do k = 1, pver
          do i = 1, ncol
             ptend%u(i,k) = utgw(i,k)
             ptend%v(i,k) = vtgw(i,k)
             ptend%s(i,k) = ttgw(i,k)
          end do
       end do


! C.-C. Chen, momentum & energy conservation
       call momentum_energy_conservation(ncol, taue, tauw, taun, taus, &
            state%pint, state%pdel, state%t, state%u, state%v, ptend%u, ptend%v, ptend%s, &
            utgw, vtgw, ttgw, src_level, dt)

! Write output fields to history file
       call outfld ('BUTGWSPEC', utgw , pcols, lchnk)
       call outfld ('BVTGWSPEC', vtgw , pcols, lchnk)
       call outfld ('BTAUE', taue , pcols, lchnk)
       call outfld ('BTAUW', tauw , pcols, lchnk)
       call outfld ('BTAUN', taun , pcols, lchnk)
       call outfld ('BTAUS', taus , pcols, lchnk)


!-Mods for Beres04

! Determine the wave source for C&M background spectrum 

       call gw_bgnd (lchnk          , ncol       ,                           &
            state%u    , state%v    , state%t    , state%pmid , state%pint , &
            state%pdel , state%rpdel, state%lnpint,kldv       , kldvmn     , &
            ksrc       , ksrcmn     , rdpldv     , tau        , ubi        , &
            ubm        , xv         , yv         , PGWV       , kbotbg     , &
            state%frontgf      ,state%frontga     )

! set phase speeds for this chunk
       call set_phase_speeds (lchnk, ncol, PGWV,  ubi , c, .false. )

! C.-C. Chen
       src_level(:) = kbotbg

! Solve for the C&M drag profile
       call gw_drag_prof (lchnk     , ncol       ,                           &
            PGWV       , kbotbg     , ktopbg     , state%u    , state%v    , &
            state%t    , state%pint , state%pdel , state%rpdel, state%lnpint,&
            rhoi       , ni         , ti         , nm         , dt         , &
            kldv       , kldvmn     , ksrc       , ksrcmn     , rdpldv     , &
            tau        , ubi        , ubm        , xv         , yv         , &
            effgw_spec , utgw       , vtgw       , tau0x      , tau0y      , &
            kvh        ,state%pmid , state%q    , state%s    , ttgw       , qtgw , .false.  ,&
            taue       , tauw       , taun       , taus       ,c, src_level)

!Add the constituent tendencies
       do m=1, pcnst

          do k = 1, pver
             do i = 1 , ncol
                ptend%q(i,k,m) = ptend%q(i,k,m) + qtgw(i,k,m) 
             end do
          end do

       end do                          ! another constituent


! add the momentum tendencies to the output tendency arrays
       do k = 1, pver
          do i = 1, ncol
             ptend%u(i,k) = ptend%u(i,k) + utgw(i,k) 
             ptend%v(i,k) = ptend%v(i,k) + vtgw(i,k) 
             ptend%s(i,k) = ptend%s(i,k) + ttgw(i,k)
          end do
       end do

! C.-C. Chen, momentum & energy conservation
       call momentum_energy_conservation(ncol, taue, tauw, taun, taus, &
            state%pint, state%pdel, state%t, state%u, state%v, ptend%u, ptend%v, ptend%s, &
            utgw, vtgw, ttgw, src_level, dt)

! Write output fields to history file
       call outfld ('UTGWSPEC', utgw , pcols, lchnk)
       call outfld ('VTGWSPEC', vtgw , pcols, lchnk)
       call outfld ('TAUE', taue , pcols, lchnk)
       call outfld ('TAUW', tauw , pcols, lchnk)
       call outfld ('TAUNET', taue+tauw , pcols, lchnk)
       call outfld ('TAUN', taun , pcols, lchnk)
       call outfld ('TAUS', taus , pcols, lchnk)

    else

! zero net tendencies if no spectrum computed
       ptend%u = 0._r8
       ptend%v = 0._r8
       ptend%s = 0._r8
       ptend%q = 0._r8

    end if
!-----------------------------------------------------------------------------
! Orographic stationary gravity wave
!-----------------------------------------------------------------------------

! Determine the orographic wave source

    call gw_oro (lchnk, ncol,                                             &
         state%u    , state%v    , state%t    , sgh        , state%pmid , &
         state%pint , state%pdel , state%zm   , nm         , pblh       , &
         kldv       , kldvmn     , ksrc       , ksrcmn     , rdpldv     , &
         tau        , ubi        , ubm        , xv         , yv         )


! set phase speeds for this chunk
    call set_phase_speeds (lchnk, ncol, 0,  ubi , c, .false.)

! C.-C. Chen
       src_level(:) = kbotoro

! Solve for the orographic drag profile
       call gw_drag_prof (lchnk     , ncol       ,                           &
            0          , kbotoro    , ktoporo    , state%u    , state%v    , &
            state%t    , state%pint , state%pdel , state%rpdel, state%lnpint,&
            rhoi       , ni         , ti         , nm         , dt         , &
            kldv       , kldvmn     , ksrc       , ksrcmn     , rdpldv     , &
            tau        , ubi        , ubm        , xv         , yv         , &
            effgw_oro  , utgw       , vtgw       , tau0x      , tau0y      , &
            kvh        ,state%pmid  , state%q    , state%s    , ttgw       , qtgw , .false.  ,&
            taue       , tauw       , taun       , taus       , c, src_level)


! Add the orographic tendencies to the spectrum tendencies
! Compute the temperature tendency from energy conservation (includes spectrum).
    do k = 1, pver
       do i = 1, ncol
          ptend%u(i,k) = ptend%u(i,k) + utgw(i,k) * landfrac(i)
          ptend%v(i,k) = ptend%v(i,k) + vtgw(i,k) * landfrac(i)
          ptend%s(i,k) = ptend%s(i,k) + ttgw(i,k)                                 &
                       -(ptend%u(i,k) * (state%u(i,k) + ptend%u(i,k)*0.5_r8*dt)   &
                       + ptend%v(i,k) * (state%v(i,k) + ptend%v(i,k)*0.5_r8*dt))
       end do
    end do

! Load TTGW with the total heating for output
    do k=1, pver
       do i=1, ncol
          ttgw(i,k) = ptend%s(i,k) 
       end do
    end do

    do m=1, pcnst

       do k = 1, pver
          do i = 1 , ncol
             ptend%q(i,k,m) = ptend%q(i,k,m) + qtgw(i,k,m) 
          end do
       end do

    end do                          ! another constituent

! Set flags for nonzero tendencies, q not yet affected by gwd
    ptend%name  = "Gravity waves drag"
    ptend%lq(:) = .true.
    ptend%ls    = .true.
    ptend%lu    = .true.
    ptend%lv    = .true.

! Write output fields to history file
    call outfld ('UTGWORO', utgw,  pcols, lchnk)
    call outfld ('VTGWORO', vtgw,  pcols, lchnk)
    call outfld ('TTGW', ttgw / cpair,  pcols, lchnk)
    call outfld ('TAUGWX',  tau0x, pcols, lchnk)
    call outfld ('TAUGWY',  tau0y, pcols, lchnk)
    call outfld ('SGH    ', sgh,   pcols, lchnk)

    return
  end  subroutine gw_intr

!===============================================================================
       subroutine momentum_energy_conservation(ncol,taue, tauw, taun, taus, &
            pi, pdel, t, u, v, dudt, dvdt, dsdt, &
            utgw, vtgw, ttgw, src_level, dt)

! C.-C. Chen, momentum & energy conservation
       implicit none

       integer, intent(in) :: ncol
       real(r8), intent(in) :: taue(pcols,0:pver)
       real(r8), intent(in) :: tauw(pcols,0:pver)
       real(r8), intent(in) :: taun(pcols,0:pver)
       real(r8), intent(in) :: taus(pcols,0:pver)
       real(r8), intent(in) :: pi(pcols,1:pver)
       real(r8), intent(in) :: pdel(pcols,pver)
       real(r8), intent(in) :: t(pcols,pver)
       real(r8), intent(in) :: u(pcols,pver)
       real(r8), intent(in) :: v(pcols,pver)
       real(r8), intent(inout) :: dudt(pcols,pver)
       real(r8), intent(inout) :: dvdt(pcols,pver)
       real(r8), intent(inout) :: dsdt(pcols,pver)
       real(r8), intent(inout) :: utgw(pcols,pver)
       real(r8), intent(inout) :: vtgw(pcols,pver)
       real(r8), intent(inout) :: ttgw(pcols,pver)
       integer, intent(in) :: src_level(pcols)
       real(r8), intent(in) :: dt
! local variables
       integer :: i,l,k
       real(r8) :: dz, dE, ut_dz, vt_dz

       do i = 1, ncol

! total mass from ground to source level, rho*dz = dp/g
        dz = 0.0_r8
        do k=pver,src_level(i)+1,-1
         dz=dz+pdel(i,k)/g
        enddo 

! tendency for U & V below source level       
        do k=src_level(i)+1,pver
         ut_dz = -(taue(i,src_level(i))+tauw(i,src_level(i)))/dz
         vt_dz = -(taun(i,src_level(i))+taus(i,src_level(i)))/dz
         dudt(i,k)=dudt(i,k)+ut_dz
         dvdt(i,k)=dvdt(i,k)+vt_dz

         utgw(i,k) = utgw(i,k)+ut_dz
         vtgw(i,k) = vtgw(i,k)+vt_dz
        enddo

! net gain/loss of total energy in the column
        dE = 0.0_r8
        do k=1,pver
         dE=dE+pdel(i,k)*(dsdt(i,k)+dudt(i,k)*(u(i,k)+dudt(i,k)*0.5_r8*dt)+  &
                                    dvdt(i,k)*(v(i,k)+dvdt(i,k)*0.5_r8*dt) )
        enddo

        dE = dE/(pi(i,pver)-pi(i,src_level(i)))

! subtract net gain/loss of total energy below source level
        do k=src_level(i)+1,pver
         dsdt(i,k) = dsdt(i,k)-dE
         ttgw(i,k) = ttgw(i,k)-dE
        enddo
       enddo
             
       return 
       end subroutine momentum_energy_conservation

!===============================================================================
  subroutine gw_prof (lchnk, ncol, u, v, t, pm, pi, rhoi, ni, ti, nm )

    use phys_grid, only: get_lat_p, get_lon_p, get_rlat_p,get_rlon_p
!-----------------------------------------------------------------------
! Compute profiles of background state quantities for the multiple
! gravity wave drag parameterization.
! 
! The parameterization is assumed to operate only where water vapor 
! concentrations are negligible in determining the density.
!-----------------------------------------------------------------------
!------------------------------Arguments--------------------------------
    integer, intent(in) :: lchnk                  ! chunk identifier
    integer, intent(in) :: ncol                   ! number of atmospheric columns

    real(r8), intent(in) :: u(pcols,pver)         ! midpoint zonal wind
    real(r8), intent(in) :: v(pcols,pver)         ! midpoint meridional wind
    real(r8), intent(in) :: t(pcols,pver)         ! midpoint temperatures
    real(r8), intent(in) :: pm(pcols,pver)        ! midpoint pressures
    real(r8), intent(in) :: pi(pcols,0:pver)      ! interface pressures

    real(r8), intent(out) :: rhoi(pcols,0:pver)   ! interface density
    real(r8), intent(out) :: ni(pcols,0:pver)     ! interface Brunt-Vaisalla frequency
    real(r8), intent(out) :: ti(pcols,0:pver)     ! interface temperature
    real(r8), intent(out) :: nm(pcols,pver)       ! midpoint brunt-vaisalla frequency

!---------------------------Local storage-------------------------------
    integer :: i,k                                ! loop indexes

    real(r8) :: dtdp
    real(r8) :: n2                                ! Brunt-Vaisalla frequency squared
    real(r8) :: ti1(pcols)                        ! Temperature at top interface
    real(r8) :: tm0(pcols)                        ! Temperature outside the model domain
    real(r8) :: pm0(pcols)                        ! pressure at which tm0 is provided (in hPa)
    real(r8) :: lat(pcols)                        ! lat indexes for current column
    real(r8) :: lon(pcols)                        ! lon indexes for current column

!-----------------------------------------------------------------------------
! Determine the interface densities and Brunt-Vaisala frequencies.
!-----------------------------------------------------------------------------

! The top interface values are calculated assuming an isothermal atmosphere 
! above the top level.
    ti1(:ncol) = t(:ncol,1)
    k = 0
    do i = 1, ncol
       ti(i,k)   = ti1(i)
       rhoi(i,k) = pi(i,k) / (r*ti(i,k))
       ni(i,k)   = sqrt (g*g / (cpair*ti(i,k)))
    end do

! Interior points use centered differences
    do k = 1, pver-1
       do i = 1, ncol
          ti(i,k)   = 0.5_r8 * (t(i,k) + t(i,k+1))
          rhoi(i,k) = pi(i,k) / (r*ti(i,k))
          dtdp      = (t(i,k+1) - t(i,k)) / (pm(i,k+1) - pm(i,k))
          n2        = g*g/ti(i,k) * (1._r8/cpair - rhoi(i,k)*dtdp)
          ni(i,k)   = sqrt (max (n2min, n2))
       end do
    end do

! Bottom interface uses bottom level temperature, density; next interface
! B-V frequency.
    k = pver
    do i = 1, ncol
       ti(i,k)   = t(i,k)
       rhoi(i,k) = pi(i,k) / (r*ti(i,k))
       ni(i,k)   = ni(i,k-1)
    end do

!-----------------------------------------------------------------------------
! Determine the midpoint Brunt-Vaisala frequencies.
!-----------------------------------------------------------------------------
    do k = 1,pver
       do i = 1,ncol
          nm(i,k) = 0.5_r8 * (ni(i,k-1) + ni(i,k))
       end do
    end do

    return
  end subroutine gw_prof

!===============================================================================

  subroutine gw_oro (lchnk, ncol,                &
       u, v, t, sgh, pm, pi, dpm, zm, nm, pblh,  &
       kldv, kldvmn, ksrc, ksrcmn, rdpldv, tau, ubi, ubm, xv, yv)
!-----------------------------------------------------------------------
! Orographic source for multiple gravity wave drag parameterization.
! 
! The stress is returned for a single wave with c=0, over orography.
! For points where the orographic variance is small (including ocean),
! the returned stress is zero. 
!------------------------------Arguments--------------------------------
    integer, intent(in) :: lchnk                  ! chunk identifier
    integer, intent(in) :: ncol                   ! number of atmospheric columns

    real(r8), intent(in) :: u(pcols,pver)         ! midpoint zonal wind
    real(r8), intent(in) :: v(pcols,pver)         ! midpoint meridional wind
    real(r8), intent(in) :: t(pcols,pver)         ! midpoint temperatures
    real(r8), intent(in) :: sgh(pcols)            ! standard deviation of orography
    real(r8), intent(in) :: pm(pcols,pver)        ! midpoint pressures
    real(r8), intent(in) :: pi(pcols,0:pver)      ! interface pressures
    real(r8), intent(in) :: dpm(pcols,pver)       ! midpoint delta p (pi(k)-pi(k-1))
    real(r8), intent(in) :: zm(pcols,pver)        ! midpoint heights
    real(r8), intent(in) :: nm(pcols,pver)        ! midpoint Brunt-Vaisalla frequency
    real(r8), intent(in) :: pblh(pcols)           ! planetary boundary layer height

    integer, intent(out) :: kldv(pcols)           ! top interface of low level stress div region
    integer, intent(out) :: kldvmn                ! min value of kldv
    integer, intent(out) :: ksrc(pcols)           ! index of top interface of source region
    integer, intent(out) :: ksrcmn                ! min value of ksrc

    real(r8), intent(out) :: rdpldv(pcols)        ! 1/dp across low level divergence region
    real(r8), intent(out) :: tau(pcols,-pgwv:pgwv,0:pver)! wave Reynolds stress
    real(r8), intent(out) :: ubi(pcols,0:pver)    ! projection of wind at interfaces
    real(r8), intent(out) :: ubm(pcols,pver)      ! projection of wind at midpoints
    real(r8), intent(out) :: xv(pcols)            ! unit vectors of source wind (x)
    real(r8), intent(out) :: yv(pcols)            ! unit vectors of source wind (y)

!---------------------------Local storage-------------------------------
    integer :: i,k                                ! loop indexes

    real(r8) :: pil                               ! don't we have pi somewhere?
    real(r8) :: lzsrc                             ! vertical wavelength at source
    real(r8) :: hdsp(pcols)                       ! surface streamline displacment height (2*sgh)
    real(r8) :: sghmax                            ! max orographic sdv to use
    real(r8) :: tauoro(pcols)                     ! c=0 stress from orography
    real(r8) :: zldv(pcols)                       ! top of the low level stress divergence region
    real(r8) :: nsrc(pcols)                       ! b-f frequency averaged over source region
    real(r8) :: psrc(pcols)                       ! interface pressure at top of source region
    real(r8) :: rsrc(pcols)                       ! density averaged over source region
    real(r8) :: usrc(pcols)                       ! u wind averaged over source region
    real(r8) :: vsrc(pcols)                       ! v wind averaged over source region

!---------------------------------------------------------------------------
! Average the basic state variables for the wave source over the depth of
! the orographic standard deviation. Here we assume that the apropiate
! values of wind, stability, etc. for determining the wave source are 
! averages over the depth of the atmosphere pentrated by the typical mountain.
! Reduces to the bottom midpoint values when sgh=0, such as over ocean.
!---------------------------------------------------------------------------
    k = pver
    do i = 1, ncol
       ksrc(i) = k-1
       psrc(i) = pi(i,k-1)
       rsrc(i) = pm(i,k)/(r*t(i,k)) * dpm(i,k)
       usrc(i) = u(i,k) * dpm(i,k)
       vsrc(i) = v(i,k) * dpm(i,k)
       nsrc(i) = nm(i,k)* dpm(i,k)
       hdsp(i) = 2.0_r8 * sgh(i)
    end do
    do k = pver-1, pver/2, -1
       do i = 1, ncol
          if (hdsp(i) > sqrt(zm(i,k)*zm(i,k+1))) then
             ksrc(i) = k-1
             psrc(i) = pi(i,k-1)
             rsrc(i) = rsrc(i) + pm(i,k) / (r*t(i,k))* dpm(i,k)
             usrc(i) = usrc(i) + u(i,k) * dpm(i,k)
             vsrc(i) = vsrc(i) + v(i,k) * dpm(i,k)
             nsrc(i) = nsrc(i) + nm(i,k)* dpm(i,k)
          end if
       end do
    end do
    do i = 1, ncol
       rsrc(i) = rsrc(i) / (pi(i,pver) - psrc(i))
       usrc(i) = usrc(i) / (pi(i,pver) - psrc(i))
       vsrc(i) = vsrc(i) / (pi(i,pver) - psrc(i))
       nsrc(i) = nsrc(i) / (pi(i,pver) - psrc(i))

       ubi(i,pver) = sqrt (usrc(i)**2 + vsrc(i)**2)
       xv(i) = usrc(i) / ubi(i,pver)
       yv(i) = vsrc(i) / ubi(i,pver)
    end do

! Project the local wind at midpoints onto the source wind.
    do k = 1, pver
       do i = 1, ncol
          ubm(i,k) = u(i,k) * xv(i) + v(i,k) * yv(i)
       end do
    end do

! Compute the interface wind projection by averaging the midpoint winds.
! Use the top level wind at the top interface.
    do i = 1, ncol
       ubi(i,0) = ubm(i,1)
    end do
    do k = 1, pver-1
       do i = 1, ncol
          ubi(i,k) = 0.5_r8 * (ubm(i,k) + ubm(i,k+1))
       end do
    end do

!---------------------------------------------------------------------------
! Determine the depth of the low level stress divergence region, as
! the max of the source region depth and 1/2 the vertical wavelength at the
! source. 
!---------------------------------------------------------------------------
    pil = acos( -1._r8 )
    do i = 1, ncol
       lzsrc   = max(2._r8 * pil * usrc(i) / nsrc(i), lzmax)
       zldv(i) = max(hdsp(i), 0.5_r8 * lzsrc)
    end do

! find the index of the top of the low level divergence region
    kldv(:) = pver-1
    do k = pver-1, pver/2, -1
       do i = 1, ncol
          if (zldv(i) .gt. sqrt(zm(i,k)*zm(i,k+1))) then
             kldv(i)  = k-1
          end if
       end do
    end do

! Determine the orographic c=0 source term following McFarlane (1987).
! Set the source top interface index to pver, if the orographic term is zero.
    do i = 1, ncol
       if ((ubi(i,pver) .gt. orovmin) .and. (hdsp(i) .gt. orohmin)) then
          sghmax = fcrit2 * (ubi(i,pver) / nsrc(i))**2
          tauoro(i) = oroko2 * min(hdsp(i)**2, sghmax) * rsrc(i) * nsrc(i) * ubi(i,pver)
       else
          tauoro(i) = 0._r8
          ksrc(i) = pver
          kldv(i) = pver
       end if
    end do

! Set the phase speeds and wave numbers in the direction of the source wind.
! Set the source stress magnitude (positive only, note that the sign of the 
! stress is the same as (c-u).
    do i = 1, ncol
       tau(i,0,pver) = tauoro(i)
    end do

! Determine the min value of kldv and ksrc for limiting later loops
! and the pressure at the top interface of the low level stress divergence
! region.

    ksrcmn = pver
    kldvmn = pver
    do i = 1, ncol
       ksrcmn = min(ksrcmn, ksrc(i))
       kldvmn = min(kldvmn, kldv(i))
       if (kldv(i) .ne. pver) then
          rdpldv(i) = 1._r8 / (pi(i,kldv(i)) - pi(i,pver))
       end if
    end do
    if (fracldv .le. 0._r8) kldvmn = pver

    return
  end subroutine gw_oro

!===============================================================================
  subroutine gw_bgnd (lchnk, ncol, &
       u, v, t, pm, pi, dpm, rdpm, piln,     &
       kldv, kldvmn, ksrc, ksrcmn, rdpldv, tau, ubi, ubm, xv, yv, &
       ngwv, kbot, frontgf, frontga)
!-----------------------------------------------------------------------
! Driver for multiple gravity wave drag parameterization.
! 
! The parameterization is assumed to operate only where water vapor 
! concentrations are negligible in determining the density.
!-----------------------------------------------------------------------

    use phys_grid,     only: get_rlat_all_p
    use physconst,     only: pi_g => pi
!   use time_manager, only: get_nstep,get_step_size,get_curr_calday	! rrg 1/7/2005

!------------------------------Arguments--------------------------------
    integer, intent(in) :: lchnk                  ! chunk identifier
    integer, intent(in) :: ncol                   ! number of atmospheric columns

    integer, intent(in) :: kbot                   ! index of bottom (source) interface
    integer, intent(in) :: ngwv                   ! number of gravity waves to use

    real(r8), intent(in) :: u(pcols,pver)         ! midpoint zonal wind
    real(r8), intent(in) :: v(pcols,pver)         ! midpoint meridional wind
    real(r8), intent(in) :: t(pcols,pver)         ! midpoint temperatures
    real(r8), intent(in) :: pm(pcols,pver)        ! midpoint pressures
    real(r8), intent(in) :: pi(pcols,0:pver)      ! interface pressures
    real(r8), intent(in) :: dpm(pcols,pver)       ! midpoint delta p (pi(k)-pi(k-1))
    real(r8), intent(in) :: rdpm(pcols,pver)      ! 1. / (pi(k)-pi(k-1))
    real(r8), intent(in) :: piln(pcols,0:pver)    ! ln(interface pressures)

    real(r8), intent(in) :: frontgf  (pcols,pver)        ! frontogenesis function
    real(r8), intent(in) :: frontga  (pcols,pver)        ! frontogenesis angle

    integer, intent(out) :: kldv(pcols)           ! top interface of low level stress divergence region
    integer, intent(out) :: kldvmn                ! min value of kldv
    integer, intent(out) :: ksrc(pcols)           ! index of top interface of source region
    integer, intent(out) :: ksrcmn                ! min value of ksrc

    real(r8), intent(in) :: rdpldv(pcols)        ! 1/dp across low level divergence region
    real(r8), intent(out) :: tau(pcols,-ngwv:ngwv,0:pver)! wave Reynolds stress
    real(r8), intent(out) :: ubi(pcols,0:pver)    ! projection of wind at interfaces
    real(r8), intent(out) :: ubm(pcols,pver)      ! projection of wind at midpoints
    real(r8), intent(out) :: xv(pcols)            ! unit vectors of source wind (x)
    real(r8), intent(out) :: yv(pcols)            ! unit vectors of source wind (y)

!---------------------------Local storage-------------------------------
    integer :: i,k,l,m                            ! loop indexes

    real(r8) :: rlat(pcols)                       ! latitude in radians for columns
    real(r8) :: tauback(pcols)                    ! background stress at c=0
    real(r8) :: flat_gw                           ! The actual lat dependence of GW spec.

    real(r8) :: usrc(pcols)                       ! u wind averaged over source region
    real(r8) :: vsrc(pcols)                       ! v wind averaged over source region

    real(r8) :: xfrontg(pcols,pver) 

!---------------------------------------------------------------------------
! Determine the source layer wind and unit vectors, then project winds.
!---------------------------------------------------------------------------

! Just use the source level interface values for the source
! wind speed and direction (unit vector).
    k = kbot
    do i = 1, ncol
       ksrc(i) = k
       kldv(i) = k
       usrc(i) = 0.5_r8*(u(i,k+1)+u(i,k))
       vsrc(i) = 0.5_r8*(v(i,k+1)+v(i,k))
       ubi(i,kbot) = sqrt (usrc(i)**2 + vsrc(i)**2)
       xv(i) = usrc(i) / ubi(i,k)
       yv(i) = vsrc(i) / ubi(i,k)
    end do

! Project the local wind at midpoints onto the source wind.
    do k = 1, kbot
       do i = 1, ncol
          ubm(i,k) = u(i,k) * xv(i) + v(i,k) * yv(i)
       end do
    end do

! Compute the interface wind projection by averaging the midpoint winds.
! Use the top level wind at the top interface.
    do i = 1, ncol
       ubi(i,0) = ubm(i,1)
    end do
    do k = 1, kbot-1
       do i = 1, ncol
          ubi(i,k) = 0.5_r8 * (ubm(i,k) + ubm(i,k+1))
       end do
    end do

   xfrontg(:,:) = 0._r8

   do i=1, ncol
      do k=1, pver
         xfrontg(i,k) = frontgf(i,k) 
      end do
   end do

   call outfld ('FRONTGF', xfrontg, pcols, lchnk)
   call outfld ('FRONTGFA', frontga, pcols, lchnk)

!-----------------------------------------------------------------------
! Gravity wave sources
!-----------------------------------------------------------------------

    tau(:,:,:) = 0._r8

! Determine the background stress at c=0
    do i=1,ncol
      tauback(i) = taubgnd * tauscal
    enddo

    do l = 1, ngwv
       do i = 1, ncol

          !
          ! GW generation depends  on frontogenesis at src level
          ! rrg: use frontgf 2 levels below launch (~600 mb) (note: hardwiring needs to be removed)
          !
          if ( xfrontg(i,kbot+2) > frontgfc) then	
             flat_gw = 1.0_r8
          else
             flat_gw = 0.0_r8
          end if
          tau(i, l,kbot) = tauback(i) * fav(l) * flat_gw
          tau(i,-l,kbot) = tau(i,l,kbot)

       end do
    end do

    do i = 1, ncol
       tau(i,0,kbot) = 0._r8
    end do

    
! Determine the min value of kldv and ksrc for limiting later loops
! and the pressure at the top interface of the low level stress divergence
! region.

    ksrcmn = pver
    kldvmn = pver

    return

  end  subroutine gw_bgnd

!===============================================================================
  subroutine gw_drag_prof (lchnk, ncol,                                  &
       ngwv,   kbot,   ktop,   u,      v,       t,                       &
       pi,     dpm,    rdpm,   piln,   rhoi,    ni,   ti,    nm,   dt,   &
       kldv,   kldvmn, ksrc,   ksrcmn, rdpldv,  tau,  ubi,   ubm,  xv,   &
       yv,     effgw,  ut,     vt,     tau0x,   tau0y,                   &
       kvh,    pm,     q,      dse,    ttgw,    qtgw, bflag ,      &
       taue,   tauw,   taun,   taus,   c, src_level)
    use diffusion_solver, only: nbot_molec
    use phys_grid,     only: get_rlat_all_p
!-----------------------------------------------------------------------
! Solve for the drag profile from the multiple gravity wave drag
! parameterization.
! 1. scan up from the wave source to determine the stress profile
! 2. scan down the stress profile to determine the tendencies
!     => apply bounds to the tendency
!          a. from wkb solution
!          b. from computational stability constraints
!     => adjust stress on interface below to reflect actual bounded tendency
!-----------------------------------------------------------------------
!------------------------------Arguments--------------------------------
    integer, intent(in) :: lchnk                  ! chunk identifier
    integer, intent(in) :: ncol                   ! number of atmospheric columns
    integer, intent(in) :: kbot                   ! index of bottom (source) interface
    integer, intent(in) :: ktop                   ! index of top interface of gwd region
    integer, intent(in) :: ngwv                   ! number of gravity waves to use
    integer, intent(in) :: kldv(pcols)            ! top interface of low level stress  divergence region
    integer, intent(in) :: kldvmn                 ! min value of kldv
    integer, intent(in) :: ksrc(pcols)            ! index of top interface of source region
    integer, intent(in) :: ksrcmn                 ! min value of ksrc

    real(r8), intent(in) :: u(pcols,pver)         ! midpoint zonal wind
    real(r8), intent(in) :: v(pcols,pver)         ! midpoint meridional wind
    real(r8), intent(in) :: t(pcols,pver)         ! midpoint temperatures
    real(r8), intent(in) :: pi(pcols,0:pver)      ! interface pressures
    real(r8), intent(in) :: dpm(pcols,pver)       ! midpoint delta p (pi(k)-pi(k-1))
    real(r8), intent(in) :: rdpm(pcols,pver)      ! 1. / (pi(k)-pi(k-1))
    real(r8), intent(in) :: piln(pcols,0:pver)    ! ln(interface pressures)
    real(r8), intent(in) :: rhoi(pcols,0:pver)    ! interface density
    real(r8), intent(in) :: ni(pcols,0:pver)      ! interface Brunt-Vaisalla frequency
    real(r8), intent(in) :: ti(pcols,0:pver)      ! interface temperature
    real(r8), intent(in) :: nm(pcols,pver)        ! midpoint Brunt-Vaisalla frequency
    real(r8), intent(in) :: dt                    ! time step
    real(r8), intent(in) :: rdpldv(pcols)         ! 1/dp across low level divergence region
    real(r8), intent(in) :: ubi(pcols,0:pver)     ! projection of wind at interfaces
    real(r8), intent(in) :: ubm(pcols,pver)       ! projection of wind at midpoints
    real(r8), intent(in) :: xv(pcols)             ! unit vectors of source wind (x)
    real(r8), intent(in) :: yv(pcols)             ! unit vectors of source wind (y)
    real(r8), intent(in) :: effgw                 ! tendency efficiency

    real(r8), intent(inout) :: tau(pcols,-pgwv:pgwv,0:pver)! wave Reynolds stress

    real(r8), intent(out) :: ut(pcols,pver)       ! zonal wind tendency
    real(r8), intent(out) :: vt(pcols,pver)       ! meridional wind tendency
    real(r8), intent(out) :: tau0x(pcols)         ! c=0 sfc. stress (zonal)
    real(r8), intent(out) :: tau0y(pcols)         ! c=0 sfc. stress (meridional)

    real(r8), intent(in) ::  kvh(pcols,0:pver)                 ! molecular thermal diffusivity
    real(r8), intent(in) ::  pm(pcols,pver)                    ! pressure at midpoints
    real(r8), intent(in) ::  q(pcols,pver,pcnst)               ! constituent array
    real(r8), intent(in) ::  dse(pcols,pver)                   ! dry static energy
    real(r8), intent(in)  :: c(pcols,-pgwv:pgwv)                       ! wave phase speeds for each column
    real(r8), intent(out) :: ttgw(pcols,pver)                  ! GW's heating
    real(r8), intent(out) :: qtgw(pcols,pver,pcnst)            ! GW's constituent tendency
    real(r8), intent(out) :: taue(pcols,0:pver)                ! Reynolds stress for eastward propagating waves
    real(r8), intent(out) :: tauw(pcols,0:pver)                ! Reynolds stress for westward propagating waves
    real(r8), intent(out) :: taun(pcols,0:pver)                ! Reynolds stress for northward propagating waves
    real(r8), intent(out) :: taus(pcols,0:pver)                ! Reynolds stress for southward propagating waves
! C.-C. Chen
    integer, intent(in) :: src_level(pcols)

!+Mods for Beres04
    logical, intent(in) :: bflag
!-Mods for Beres04

!---------------------------local storage-------------------------------
    integer :: i,k,l,m,n                          ! loop indexes

    real(r8) :: d(pcols)                          ! "total" diffusivity 
    real(r8) :: dsat(pcols,-ngwv:ngwv)            ! saturation diffusivity
    real(r8) :: dscal                             ! fraction of dsat to use
    real(r8) :: mi                                ! imaginary part of vertical wavenumber
    real(r8) :: taudmp                            ! stress after damping
    real(r8) :: taumax(pcols)                     ! max(tau) for any l
    real(r8) :: tausat(pcols,-ngwv:ngwv)          ! saturation stress
    real(r8) :: ubmc(pcols,-ngwv:ngwv)            ! (ub-c)
    real(r8) :: ubmc2                             ! (ub-c)**2
    real(r8) :: ubt(pcols,pver)                   ! ubar tendency
    real(r8) :: ubtl                              ! ubar tendency from wave l
    real(r8) :: ubtlsat                           ! saturation tendency
    real(r8) :: wrk
    real(r8) :: gwut(pcols,pver,-ngwv:ngwv)       ! gravity wave wind tendency for each wave
    real(r8) :: egwdffi(pcols,0:pver)             ! effective gw diffusivity at interfaces

! C.-C. Chen
    real(r8) :: tausg(pcols,-ngwv:ngwv,0:pver)            ! signed wave Reynolds stress 

    real(r8) :: taulu(pcols,0:pver)                      ! Reynolds stress for waves propagating to the left
                                                  ! of the source zonal winds
    real(r8) :: tauru(pcols,0:pver)                      ! Reynolds stress for waves propagating to the right
                                                  ! of the source zonal winds
    real(r8) :: taulv(pcols,0:pver)                      ! Reynolds stress for waves propagating to the left
                                                  ! of the source meridional winds
    real(r8) :: taurv(pcols,0:pver)                      ! Reynolds stress for waves propagating to the right
                                                  ! of the source meridional winds
    real(r8) :: ut1(pcols,pver)                   ! zonal wind tendency for  c < -40
    real(r8) :: ut2(pcols,pver)                   ! zonal wind tendency for  -40 < c < -15
    real(r8) :: ut3(pcols,pver)                   ! zonal wind tendency for  -15 < c <  15
    real(r8) :: ut4(pcols,pver)                   ! zonal wind tendency for   15 < c <  40
    real(r8) :: ut5(pcols,pver)                   ! zonal wind tendency for   40 < c


    real(r8) :: EMF(pcols,pver)
    real(r8) :: WMF(pcols,pver)
    real(r8) :: NMF(pcols,pver)
    real(r8) :: SMF(pcols,pver)

    real(r8) :: cx,cy                             ! project wave phase speeds

    real(r8) :: rlat(pcols)                       ! latitude of the columns in the chunck (radians)

    real(r8) :: ptaper(pcols)                     ! rrg: polar taper
    real(r8) :: al0				  ! rrg: used in lat dependence of the taper
    real(r8) :: dlat0                             ! rrg: used in lat dependence of the taper


    real(r8) :: dummy(pcols,pver)
    real(r8) :: dummyx(pcols,pver)
    real(r8) :: dummyy(pcols,pver)
    real(r8) :: ccn
    character(len=80) fnum,fmt ,dumc1, dumc1x, dumc1y          ! dummies

    real(r8):: taux(pcols,-pgwv:pgwv,0:pver)      ! wave stress in zonal direction
    real(r8):: tauy(pcols,-pgwv:pgwv,0:pver)      ! wave stress in meridional direction

    integer  :: ix, iy                            ! dummies
    integer  :: usi,vsi                           ! dummies
    character(len=80) :: dum1,dum2,dum3,dum4      ! dummies
 
! Initialize gravity wave drag tendencies to zero

!+Mods for Beres04
    ut1(:,:) = 0._r8
    ut2(:,:) = 0._r8
    ut3(:,:) = 0._r8
    ut4(:,:) = 0._r8
    ut5(:,:) = 0._r8
 



!-Mods for Beres04

    call get_rlat_all_p(lchnk, ncol, rlat)

    do k = 1,pver
       do i = 1,pcols
          ut(i,k) = 0._r8
          vt(i,k) = 0._r8
       end do
    end do

    taue(:,:) = 0._r8
    tauw(:,:) = 0._r8
    taun(:,:) = 0._r8
    taus(:,:) = 0._r8

    EMF(:,:)=0._r8
    WMF(:,:)=0._r8
    SMF(:,:)=0._r8
    NMF(:,:)=0._r8


    gwut(:,:,:) = 0._r8

    taux(:,:,:)=0._r8
    tauy(:,:,:)=0._r8

    if ( ngwv > 0) then

       ! Load bottom values of Reynolds stresses
!       k=kbot

       do l=-ngwv, ngwv
          do i=1, ncol
             k = src_level(i)
             tausg(i,l,k) = sign( tau(i,l,k), c(i,l)-ubi(i,k) )
             end do
          end do

          taulu(:,:) = 0._r8
          tauru(:,:) = 0._r8
          do i=1, ncol
             do l=-ngwv, ngwv

                k = src_level(i)

                if ( c(i,l) > ubi(i,k) ) then

                   tauru(i,k) = tauru(i,k) + tausg(i,l,k)

                else if ( c(i,l) < ubi(i,k) ) then

                   taulu(i,k) = taulu(i,k) + tausg(i,l,k)

                end if

             end do
          end do

          do i=1, ncol

             k = src_level(i)

             if (xv(i) > 0._r8) then

                taue(i,k) = tauru(i,k) * xv(i)
                tauw(i,k) = taulu(i,k) * xv(i)

             else if ( xv(i) < 0.) then

                taue(i,k) = taulu(i,k) * xv(i)
                tauw(i,k) = tauru(i,k) * xv(i)

             end if

             if ( yv(i) > 0._r8) then

                taun(i,k) = tauru(i,k) * yv(i)
                taus(i,k) = taulu(i,k) * yv(i)

             else if ( yv(i) < 0.) then

                taun(i,k) = taulu(i,k) * yv(i)
                taus(i,k) = tauru(i,k) * yv(i)

             end if

          end do

       end if

!---------------------------------------------------------------------------
! Compute the stress profiles and diffusivities
!---------------------------------------------------------------------------

! Loop from bottom to top to get stress profiles      

! C.-C. Chen, move k loop into i loop
!    do k = kbot-1, ktop, -1

! Determine the absolute value of the saturation stress and the diffusivity
! for each wave.
! Define critical levels where the sign of (u-c) changes between interfaces.

       do l = -ngwv, ngwv
          do i = 1, ncol

          do k = src_level(i)-1, ktop, -1

             ubmc(i,l) = ubi(i,k) - c(i,l)
             tausat(i,l) = abs (effkwv * rhoi(i,k) * ubmc(i,l)**3 / (2._r8*ni(i,k)) )
             if (tausat(i,l) .le. taumin) tausat(i,l) = 0.0_r8

             if ( ubmc(i,l) / (ubi(i,k+1) - c(i,l) + 1e-8_r8) .le. 0.0_r8) tausat(i,l) = 0.0_r8

! Compute stress for each wave. The stress at this level is the min of 
! the saturation stress and the stress at the level below reduced by damping.
! The sign of the stress must be the same as at the level below.

       if (k <= nbot_molec) then

! Total diffusivity
!          d(:ncol) = dback
!          do i=1, ncol
             d(i) = dback + kvh(i,k) 
!          end do

!          do l = -ngwv, ngwv
!             do i = 1, ncol
                ubmc2 = max(ubmc(i,l)**2, ubmc2mn)
                mi = ni(i,k) / (2._r8 * kwv * ubmc2) * (alpha(k) + ni(i,k)**2/ubmc2 * d(i))
                wrk = -2._r8*mi*rog*t(i,k+1)*(piln(i,k+1) - piln(i,k))
                if( wrk >= -150._r8 ) then
                   taudmp = tau(i,l,k+1) * exp( wrk )
                else
                   taudmp = 0._r8
                end if
                if( taudmp <= taumin ) then
		      taudmp = 0._r8
                end if
                tau(i,l,k) = min (taudmp, tausat(i,l))			
!             end do
!          end do
       else
!          do l= -ngwv, ngwv
!             do i=1, ncol
                tau(i,l,k) = min (tau(i,l,k+1), tausat(i,l))
!             end do
!          end do

       end if

! The orographic stress term does not change across the source region
! Note that k ge ksrcmn cannot occur without an orographic source term

       if((l.eq.0).and.(k.ge.ksrcmn).and.(k.ge.ksrc(i))) then
           tau(i,0,k) = tau(i,0,pver)
       endif

 !      if (k .ge. ksrcmn) then
 !         do i = 1, ncol
 !            if (k .ge. ksrc(i)) then
 !               tau(i,0,k) = tau(i,0,pver) 
 !            end if
 !         end do
 !      end if

! Require that the orographic term decrease linearly (with pressure) 
! within the low level stress divergence region. This supersedes the above
! requirment of constant stress within the source region.
! Note that k ge kldvmn cannot occur without an orographic source term, since
! kldvmn=pver then and k<=pver-1

       if((l.eq.0).and.(k.ge.kldvmn).and.(k.ge.kldv(i))) then
                tau(i,0,k) = min (tau(i,0,k), tau(i,0,pver)  * &
                     (1._r8 - fracldv * (pi(i,k) - pi(i,pver)) * rdpldv(i)))
       endif

!       if (k .ge. kldvmn) then
!          do i = 1, ncol
!             if (k .ge. kldv(i)) then
!                tau(i,0,k) = min (tau(i,0,k), tau(i,0,pver)  * &
!                     (1._r8 - fracldv * (pi(i,k) - pi(i,pver)) * rdpldv(i)))
!             end if
!          end do
!       end if

          enddo
          end do
       end do





! Diagnostic output of TAU
       if ( ngwv > 0) then
          do l=-ngwv, ngwv
             do i=1, ncol
             do k=src_level(i)-1,ktop,-1
                tausg(i,l,k) = sign( tau(i,l,k), c(i,l)-ubi(i,k) )

                taulu(i,k) = 0._r8
                tauru(i,k) = 0._r8
                taulv(i,k) = 0._r8
                taurv(i,k) = 0._r8

             enddo
             enddo
           enddo

          do l = -ngwv,ngwv
             do i = 1,ncol
             do k = src_level(i)-1,ktop,-1

                if ( c(i,l) < ubi(i,src_level(i)) ) then
                   taulu(i,k) = taulu(i,k) + tausg(i,l,k)
                   taulv(i,k) = taulv(i,k) + tausg(i,l,k)
                else if (  c(i,l) > ubi(i,src_level(i))  ) then
                   tauru(i,k) = tauru(i,k) + tausg(i,l,k)
                   taurv(i,k) = taurv(i,k) + tausg(i,l,k)
                end if
!             end do
!          end do
             
             enddo
             end do
          end do



          do i = 1, ncol
          do k = src_level(i)-1,ktop,-1
             if (xv(i) > 0._r8) then
                taue(i,k) = tauru(i,k) * xv(i)
                tauw(i,k) = taulu(i,k) * xv(i)
             else if ( xv(i) < 0._r8) then
                taue(i,k) = taulu(i,k) * xv(i)
                tauw(i,k) = tauru(i,k) * xv(i)
             end if

             if ( yv(i) > 0._r8) then
                taun(i,k) = taurv(i,k) * yv(i)
                taus(i,k) = taulv(i,k) * yv(i)
             else if ( yv(i) < 0._r8) then
                taun(i,k) = taulv(i,k) * yv(i)
                taus(i,k) = taurv(i,k) * yv(i)
             end if
          enddo
          end do
       end if

!    end do


!---------------------------------------------------------------------------
! Compute the tendencies from the stress divergence.
!---------------------------------------------------------------------------

!++rrg: define polar taper 

  al0   = 82.5_r8 * d2r ! 3.14159_r8/180._r8
  dlat0 =  5.0_r8 * d2r ! 3.14159_r8/180._r8

  if(ngwv > 0 .and. .not.bflag) then 	! taper CM only
    do i = 1, ncol
!     ptaper(i) = 0.25_r8 * (1.+tanh((rlat(i)+al0)/dlat0)) * (1.-tanh((rlat(i)-al0)/dlat0)) 
      ptaper(i) = cos(rlat(i))
    end do
  else					! do not taper other cases
      ptaper(:) = 1._r8
  endif

!--rrg

!----------------------------------------------------------------------------
!  do k = ktop+1, kbot			! Loop over levels from top to bottom
!----------------------------------------------------------------------------

! Accumulate the mean wind tendency over wavenumber.
       do i = 1, ncol
       do k = ktop+1,src_level(i)
          ubt(i,k) = 0.0_r8
       enddo
       end do

       do l = -ngwv, ngwv	! loop over wave
          do i = 1, ncol	! loop over column
          do k = ktop+1,src_level(i)

! Determine the wind tendency including excess stress carried down from above.
             ubtl = g * (tau(i,l,k)-tau(i,l,k-1)) * rdpm(i,k)

! Apply tendency limits to maintain numerical stability.
! 1. du/dt < |c-u|/dt  so u-c cannot change sign (u^n+1 = u^n + du/dt * dt)
! 2. du/dt < tndmax    so that ridicuously large tendencies are not permitted
             ubtl = min(ubtl, umcfac * abs(c(i,l)-ubm(i,k)) / dt)
             ubtl = min(ubtl, tndmax) 

! accumulate the mean wind tendency over wavenumber, aplying efficiency and taper:
             ubt (i,k) = ubt (i,k) + sign(ubtl, c(i,l)-ubm(i,k)) * effgw * ptaper(i)

! save tendency for each wave (for later computation of kzz), applying efficiency and taper:
             gwut(i,k,l) = sign(ubtl, c(i,l)-ubm(i,k)) * effgw * ptaper(i)

! Accumulate wind tendencies based on the phase speed of spectral components
             if (ngwv > 0 ) then

                IF (c(i,l).lt.-40._r8) THEN

                   ut1(i,k) = ut1(i,k) + sign(ubtl, c(i,l)-ubm(i,k)) * effgw * ptaper(i)

                ELSE IF ((c(i,l).ge.-40._r8).and.(c(i,l).lt.-15._r8)) THEN

                   ut2(i,k) = ut2(i,k) + sign(ubtl, c(i,l)-ubm(i,k)) * effgw * ptaper(i)

                ELSE IF ((c(i,l).ge.-15._r8).and.(c(i,l).le.15._r8)) THEN

                   ut3(i,k) = ut3(i,k) + sign(ubtl, c(i,l)-ubm(i,k)) * effgw * ptaper(i)

                ELSE IF ((c(i,l).gt.15._r8).and.(c(i,l).le.40._r8)) THEN

                   ut4(i,k) = ut4(i,k) + sign(ubtl, c(i,l)-ubm(i,k)) * effgw * ptaper(i)

                ELSE IF (c(i,l).gt.40._r8) THEN

                   ut5(i,k) = ut5(i,k) + sign(ubtl, c(i,l)-ubm(i,k)) * effgw * ptaper(i)

                END IF

             end if

! Redetermine the effective stress on the interface below from the wind 
! tendency. If the wind tendency was limited above, then the new stress
! will be smaller than the old stress and will cause stress divergence in
! the next layer down. This has the effect of smoothing large stress 
! divergences downward while conserving total stress.

             tau(i,l,k) = tau(i,l,k-1) + ubtl * dpm(i,k) / g

          enddo
          end do		! end loop over column
       end do			! end loop over wave


! rrg: project the mean wind tendency (modified by effgw) onto x,y components. 
       do i = 1, ncol
       do k = ktop+1, src_level(i)

           ut(i,k) = ubt(i,k) * xv(i) 
           vt(i,k) = ubt(i,k) * yv(i)

	if (ngwv > 0) then ! zonal tendencies for spectral bands

           ut1(i,k) = ut1(i,k) * xv(i) 
           ut2(i,k) = ut2(i,k) * xv(i) 
           ut3(i,k) = ut3(i,k) * xv(i)
           ut4(i,k) = ut4(i,k) * xv(i)
           ut5(i,k) = ut5(i,k) * xv(i)


        endif

       enddo
       end do

!---------------------------------------------------------------------------
!  end do 			! end loop over level 
!---------------------------------------------------------------------------


    if (ngwv > 0) then

       if (bflag) then

          ! Output wave tendencies binned by phase speed

          call outfld ('BUTEND1', ut1, pcols, lchnk)
          call outfld ('BUTEND2', ut2, pcols, lchnk)
          call outfld ('BUTEND3', ut3, pcols, lchnk)
          call outfld ('BUTEND4', ut4, pcols, lchnk)
          call outfld ('BUTEND5', ut5, pcols, lchnk)


       else

          call outfld ('UTEND1', ut1, pcols, lchnk)
          call outfld ('UTEND2', ut2, pcols, lchnk)
          call outfld ('UTEND3', ut3, pcols, lchnk)
          call outfld ('UTEND4', ut4, pcols, lchnk)
          call outfld ('UTEND5', ut5, pcols, lchnk)


       end if

       call gw_ediff (lchnk   ,ncol                               &
            ,ngwv     ,ktopbg   ,kbotbg  ,gwut      ,ubm          &
            ,nm       ,t       ,ksrc     ,egwdffi  ,c) 
       
       do m = 1, pcnst

          call gw_const_tend  (lchnk   ,ncol                      &
               ,ngwv         ,kbotbg      ,ktopbg     ,q(1,1,m)   &
               ,pi          ,pm           ,t          ,dt         &
               ,rdpm        ,qtgw(1,1,m)  ,egwdffi     )

       enddo

       call gw_temp_tend  (lchnk           ,ncol                 &
            ,ngwv     ,kbotbg     ,ktopbg  ,dse                  &
            ,pi       ,pm         ,t       ,gwut                 &
            ,dt       ,rdpm       ,egwdffi ,ttgw     ,c)

    else 

       qtgw(:,:,:) = 0._r8
       ttgw(:,:) = 0._r8
       
    end if

!-----------------------------------------------------------------------
! Project the c=0 stress (scaled) in the direction of the source wind for recording
! on the output file.
!-----------------------------------------------------------------------
    do i = 1, ncol
       tau0x(i) = tau(i,0,kbot) * xv(i) * effgw
       tau0y(i) = tau(i,0,kbot) * yv(i) * effgw
    end do


! Output Tau with height - output taux and tauy separately

    if (ngwv > 0) then

       do i=1,ncol
          do l=-ngwv,ngwv
             cx = c(i,l)*xv(i)
             cy = c(i,l)*yv(i)
             ix = NINT (cx/dc)
             iy = NINT (cy/dc)
             if (abs(ix)>ngwv) then
                ix = sign(ngwv,ix)
             endif
             if (abs(iy)>ngwv) then
                iy = sign(ngwv,iy)
             endif
             taux(i,ix,1:pver) = taux(i,ix,1:pver) + ABS(tau(i,l,1:pver)*xv(i))
             tauy(i,iy,1:pver) = tauy(i,iy,1:pver) + ABS(tau(i,l,1:pver)*yv(i))
          enddo
       enddo

       do i=1,ncol
          do k=1,pver
             usi = NINT (u(i,k)/dc)  
             if (abs(usi).gt.ngwv) then
                usi=sign(ngwv,usi)
             endif

             vsi = NINT (v(i,k)/dc)
             if (abs(vsi).gt.ngwv) then
                vsi=sign(ngwv,vsi)
             endif

             do l=-ngwv,usi-1
                WMF(i,k) = WMF(i,k) + (taux(i,l,k)*effgw)
             enddo

             do l=usi,ngwv
                EMF(i,k) = EMF(i,k) + (taux(i,l,k)*effgw)
             enddo

             do l=-ngwv,vsi-1
                SMF(i,k) = SMF(i,k) + (tauy(i,l,k)*effgw)
             enddo

             do l=vsi,ngwv
                NMF(i,k) = NMF(i,k) + (tauy(i,l,k)*effgw)
             enddo
          enddo
       enddo
       
       if (bflag) then
          call outfld ('BEMF',EMF(1:pcols,1:pver),pcols,lchnk)
          call outfld ('BWMF',WMF(1:pcols,1:pver),pcols,lchnk)
          call outfld ('BNMF',NMF(1:pcols,1:pver),pcols,lchnk)
          call outfld ('BSMF',SMF(1:pcols,1:pver),pcols,lchnk)
       else
          call outfld ('EMF',EMF(1:pcols,1:pver),pcols,lchnk)
          call outfld ('WMF',WMF(1:pcols,1:pver),pcols,lchnk)
          call outfld ('NMF',NMF(1:pcols,1:pver),pcols,lchnk)
          call outfld ('SMF',SMF(1:pcols,1:pver),pcols,lchnk)
       endif

! Positive phase velocities

       if (bflag) then
          do m=0,ngwv

             dummyx(:pcols,:pver)=0._r8
             dummyy(:pcols,:pver)=0._r8

             do i=1,ncol
                do l=-ngwv,ngwv
                   if(abs(cref(m)-c(i,l)) <= dc/2.) then

                      dummyx(i,1:pver)=taux(i,l,1:pver)
                      dummyy(i,1:pver)=tauy(i,l,1:pver)

                   end if

                enddo
             enddo
             write (fnum,fmt='(i3)') 100+m

             dumc1x='BTAUXSp'//fnum(2:3)
             dumc1y='BTAUYSp'//fnum(2:3)

             call outfld (dumc1x,dummyx(1:pcols,1:pver),pcols,lchnk)
             call outfld (dumc1y,dummyy(1:pcols,1:pver),pcols,lchnk)

          enddo
       else
          do m=0,ngwv

             dummyx(:pcols,:pver)=0._r8
             dummyy(:pcols,:pver)=0._r8

             do i=1,ncol
                do l=-ngwv,ngwv
                   if(abs(cref2(m)-c(i,l)) <= dc/2.) then

                      dummyx(i,1:pver)=taux(i,l,1:pver)
                      dummyy(i,1:pver)=tauy(i,l,1:pver)

                   end if

                enddo
             enddo
             write (fnum,fmt='(i3)') 100+m

             dumc1x='TAUXSp'//fnum(2:3)
             dumc1y='TAUYSp'//fnum(2:3)

             call outfld (dumc1x,dummyx(1:pcols,1:pver),pcols,lchnk)
             call outfld (dumc1y,dummyy(1:pcols,1:pver),pcols,lchnk)

          enddo
       endif

! Negative phase velocities

       if (bflag) then	
	  n=0
          do m=-1,-ngwv, -1
             n=n+1
             dummyx(:pcols,:pver)=0._r8
             dummyy(:pcols,:pver)=0._r8

             do i=1,ncol
                do l=-ngwv,ngwv
                   if(abs(cref(m)-c(i,l)) <= dc/2.) then

                      dummyx(i,1:pver)=taux(i,l,1:pver)
                      dummyy(i,1:pver)=tauy(i,l,1:pver)


                   end if

                enddo
             enddo
             write (fnum,fmt='(i3)') 100+n

             dumc1x='BTAUXSn'//fnum(2:3)
             dumc1y='BTAUYSn'//fnum(2:3)

             call outfld (dumc1x,dummyx(1:pcols,1:pver),pcols,lchnk)
             call outfld (dumc1y,dummyy(1:pcols,1:pver),pcols,lchnk)
          enddo
       else      
	  n=0
          do m=-1,-ngwv, -1
             n=n+1
             dummyx(:pcols,:pver)=0._r8
             dummyy(:pcols,:pver)=0._r8

             do i=1,ncol
                do l=-ngwv,ngwv
                   if(abs(cref2(m)-c(i,l)) <= dc/2.) then

                      dummyx(i,1:pver)=taux(i,l,1:pver)
                      dummyy(i,1:pver)=tauy(i,l,1:pver)

                   end if

                enddo
             enddo
             write (fnum,fmt='(i3)') 100+n

             dumc1x='TAUXSn'//fnum(2:3)
             dumc1y='TAUYSn'//fnum(2:3)

             call outfld (dumc1x,dummyx(1:pcols,1:pver),pcols,lchnk)
             call outfld (dumc1y,dummyy(1:pcols,1:pver),pcols,lchnk)
          enddo
       endif
    end if  !ngwv > 0

    return
  end subroutine gw_drag_prof

!=======================================================================================
subroutine gw_const_tend (lchnk   ,ncol                                            &
                           ,ngwv    ,kbot      ,ktop    ,q       ,pi      ,pm      &
                           ,t       ,dt        ,rdpm    ,dq      ,egwdffi )

!
!  calculates constituent tendencies due to gw breaking
!
!  method:
!  a constituent flux on interfaces is given by
!
!                                   
!              rho * (w'q') = rho * Deff qz
!
!
! where (all evaluated on interfaces):
!
!        rho   = density
!        qz    = constituent vertical gradient
!        Deff  = effective diffusivity
!
!
! an effective diffusivity is calculated by adding up the diffusivities from all waves
! (see gw_ediff)
! tendency is calculated by invoking  lu decomposition and solving as for a regular
! diffusion equation.
!
! author : sassi - jan 2001
!---------------------------------------------------------------------------------

 implicit none

!---------------------------input arguments---------------------------------------
    integer, intent(in) :: lchnk                  ! chunk identifier
    integer, intent(in) :: ncol                   ! number of atmospheric columns
    integer, intent(in) :: kbot                   ! index of bottom (source) midpoint
    integer, intent(in) :: ktop                   ! index of top midpoint of gwd region
    integer, intent(in) :: ngwv                   ! number of gravity waves to use

    real(r8), intent(in) :: q(pcols,pver)                 ! constituent
    real(r8), intent(in) :: t(pcols,pver)                 ! temperature at midpoints
    real(r8), intent(in) :: dt                            ! time step
    real(r8), intent(in) :: rdpm(pcols,pver)              ! 1./pdel  (thickness bet interfaces)
    real(r8), intent(in) :: pi(pcols,pverp)               ! interface pressure
    real(r8), intent(in) :: pm(pcols,pver)                ! midpoint  pressure
    real(r8), intent(in) :: egwdffi(pcols,pverp)          ! effective gw diffusivity at interfaces

!--------------------------output arguments-------------------------------------
 
    real(r8), intent(out) :: dq(pcols,pver)              ! constituent tendencies at midpoints

!--------------------------local worspace---------------------------------------

    integer i,k,l
    real(r8) tmpi2(pcols,pverp)                  ! factor
    real(r8) ca(pcols,pver)                       ! -upper diagonal
    real(r8) cc(pcols,pver)                       ! -lower diagonal
    real(r8) term(pcols,pver)                     ! 1./(1. + ca(k) + cc(k) - cc(k)*ze(k-1))
    real(r8) ze(pcols,pver)                       ! term in tri-diag. matrix system
    real(r8) qconst(pcols,pver)                   ! temporary storage for constituent
    real(r8) rho(pcols,pverp)
    real(r8) zero(pcols)

!-------------------------------------------------------------------------------
    zero       = 0.0_r8
    dq(:,:)    = 0.0_r8

! Compute factor for diffusion matrices at interfaces (except top), 
! Compute rho at interfaces p/RT,  Ti = (Tm_k + Tm_k-1)/2,  interface k-
    do k = 2, pver
       do i = 1, ncol
          rho(i,k) = pi(i,k) * 2._r8 / (r * (t(i,k) + t(i,k-1)))
       end do
    end do
    do i = 1, ncol
       rho(i,pverp) = pi (i,pverp) / (r * t(i,pver))
    end do
    do i = 1, ncol
       rho(i,1) = pi (i,1) / (r * t(i,1))
    end do

! Calculate dt * (g*rho)^2/dp at interior interfaces,  interface k-
    tmpi2(:,:) = 0.0_r8
    do k = ktop+2 , kbot+1
       do i = 1, ncol
          tmpi2(i,k) = dt * (g * rho(i,k))**2 / (pm(i,k) - pm(i,k-1))
       end do
    end do

   qconst = 0.0_r8
   do k = 1, pver
      do i = 1, ncol
         qconst(i,k) = q(i,k)
      enddo
   enddo

! Decompose and solve the diffusion matrix
     call vd_lu_decomp(                                                &
          lchnk      ,ncol       ,                                     &
          zero       ,egwdffi    ,tmpi2      ,rdpm       ,dt         ,zero , &
          ca         ,cc         ,term       ,ze         ,ktop+1     , &
          kbot+1 )
       call vd_lu_solve(                                                 &
            lchnk      ,ncol       ,                                     &
            qconst(1,1),ca         ,ze         ,term       ,             &
            ktop+1      ,kbot+1    ,zero   )

! Evaluate tendency to be reported back
      do k = ktop+1, kbot
         do i = 1, ncol
            dq(i,k) = (qconst(i,k) - q(i,k)) / dt
         enddo
     enddo
 
  return
end subroutine gw_const_tend


!=====================================================================================

subroutine gw_ediff (lchnk   ,ncol                             &
                    ,ngwv    ,ktop    ,kbot   ,gwut    ,ubm    &
                    ,nm      ,t                       &
                    ,ksrc    ,egwdffi ,c        )

! 
! Calculate effective diffusivity associated with GW forcing.
!
! Author: F. Sassi, Jan 31, 2001
!

implicit none

!-------------------------------Input Arguments------------------------------------------

    integer, intent(in) :: lchnk                  ! chunk identifier
    integer, intent(in) :: ncol                   ! number of atmospheric columns
    integer, intent(in) :: kbot                   ! index of bottom (source) midpoint
    integer, intent(in) :: ktop                   ! index of top midpoint of gwd region
    integer, intent(in) :: ngwv                   ! number of gravity waves to use
    integer, intent(in) :: ksrc(pcols)            ! index of top interface of source region

    real(r8), intent(in) :: t(pcols,pver)         ! midpoint temperature
    real(r8), intent(in) :: gwut(pcols,pver,-ngwv:ngwv)   ! GW zonal wind tendencies at midpoint
    real(r8), intent(in) :: ubm(pcols,pver)               ! projection of wind at midpoints
    real(r8), intent(in) :: nm(pcols,pver)                ! Brunt-Vaisalla frequency    
    real(r8), intent(in) :: c(pcols,-ngwv:ngwv)           ! wave phase speeds for each column

!-----------------------------Output Argument-------------------------------------------
    real(r8), intent(out) :: egwdffi(pcols,0:pver)  ! effective gw diffusivity at interfaces

!-----------------------------Local workspace-------------------------------------------

    real(r8) egwdffm(pcols,pver)                   ! effective gw diffusivity at midpoints
    integer i,k,l,m
!   real(r8), parameter :: prndl=0.50              ! Inverse Prandtl no.       !fs
    real(r8), parameter :: prndl=0.25_r8           ! Inverse Prandtl no.       !rrg 1/25/2005
    real(r8) :: coef
    real(r8) :: nms2(pver)                         ! zonal average nm^2
    real(r8) :: rho(pcols,pverp)                   ! density at interfaces
    real(r8), parameter :: dscale=7000._r8         ! density scale height

    real(r8) :: dummy(pcols,pver)
    real(r8) :: ccn,cmx,cmn
    character(len=8) :: fnum,fname,dumc1           ! dummy characters

!---------------------------------------------------------------------------------------

! zero
  
  egwdffi(:,:) = 0._r8
  egwdffm(:,:) = 0._r8

  if (ngwv == 0) return
  

! Calculate effective diffusivity at midpoints
    do k = ktop+1, kbot
       do i = 1, ncol
          do l = -ngwv, ngwv

             egwdffm(i,k) = egwdffm(i,k) + &
                    prndl * 0.5_r8 * gwut(i,k,l) * (c(i,l)-ubm(i,k)) / nm(i,k)**2  !rrg/bab

          enddo                 ! another wave
      enddo                     ! another long
    enddo                       ! another level


! Interpolate effective diffusivity to interfaces
   do i = 1, ncol
! Assume zero  at top and bottom interfaces
      do k=ktop+1, kbot-1
         egwdffi(i,k) = .5_r8*(egwdffm(i,k) + egwdffm(i,k+1))
      enddo
   enddo

! Limit diffusivity to some reasonable value
    do i=1, ncol
       do k=0, pver
          egwdffi(i,k) = min( 150._r8,egwdffi(i,k) )
       enddo
    enddo


! Do not calculate diffusivities below source level
   do k= ktop,kbot
      do i = 1, ncol
         if (k .ge. ksrc(i)) egwdffi(i,k) = 0.0_r8
      enddo
   enddo

   fname='EKGWSPEC'
   call outfld (fname,egwdffi(1:pcols,0:pver),pcols,lchnk)

   return
end subroutine gw_ediff

!========================================================================================

  subroutine gw_temp_tend (lchnk   ,ncol                                              &
                          ,ngwv    ,kbot   ,ktop    ,dse     ,pi       ,pm            &
                          ,t       ,gwut   ,dt      ,rdpm                    &
                          ,egwdffi ,ttgw   ,c)

!
!  Calculates temperature tendency due to GW breaking
!
!  Method:
!  GW mixing is calculated using the effective diffusivity obtained from gw_const_tend: 
!  it is used to diffuse dry static energy.
!
! Author : Sassi - Jan 2001
!---------------------------------------------------------------------------------

 implicit none

!---------------------------Input Arguments---------------------------------------
    integer, intent(in) :: lchnk                  ! chunk identifier
    integer, intent(in) :: ncol                   ! number of atmospheric columns
    integer, intent(in) :: kbot                   ! index of bottom (source) interface
    integer, intent(in) :: ktop                   ! index of top interface
    integer, intent(in) :: ngwv                   ! number of gravity waves to use

    real(r8), intent(in) :: gwut(pcols,pver,-ngwv:ngwv)   ! GW zonal wind tendencies at midpoint
    real(r8), intent(in) :: t(pcols,pver)                 ! temperature at midpoints
    real(r8), intent(in) :: dt                            ! time step
    real(r8), intent(in) :: rdpm(pcols,pver)              ! 1./pdel  (thickness bet interfaces)
    real(r8), intent(in) :: pi(pcols,pverp)               ! interface pressure
    real(r8), intent(in) :: pm(pcols,pver)                ! midpoint pressure
    real(r8), intent(in) :: egwdffi(pcols,pverp)          ! effective gw diffusivity at interfaces
    real(r8), intent(in) :: dse(pcols,pver)               ! midpoint dry static energy
    real(r8), intent(in) :: c(pcols,-ngwv:ngwv)                       ! wave phase speeds for each column

!--------------------------Output Arguments-------------------------------------

    real(r8), intent(out) :: ttgw(pcols,pver) 
 
!--------------------------Local Worspace---------------------------------------

    integer i,k,m,l
    real(r8) :: delta                        ! dissipation rate
    real(r8) :: tmpi2(pcols,pverp)           ! pi(k)/(.5*(tm1(k)+tm1(k-1))
    real(r8) :: ca(pcols,pver)               ! -upper diagonal
    real(r8) :: cc(pcols,pver)               ! -lower diagonal
    real(r8) :: term(pcols,pver)             ! 1./(1. + ca(k) + cc(k) - cc(k)*ze(k-1))
    real(r8) :: ze(pcols,pver)               ! term in tri-diag. matrix system
    real(r8) :: temp(pcols,pver)             ! temporary storage for constituents
    real(r8) :: rho(pcols,pverp)
    real(r8) :: zero(pcols)
    real(r8) :: dttdf(pcols,pver)            ! dse tendencies at midpoints: diffusion term
    real(r8) :: dttke(pcols,pver)            ! dse tendencies at midpoints: kin. energy conversion term

!-------------------------------------------------------------------------------
    zero(:)    = 0.0_r8
    dttdf(:,:) = 0.0_r8
    dttke(:,:) = 0.0_r8

    if (ngwv .eq. 0) return

!
! Compute factor for diffusion matrices at interior interfaces 
! Note that [ktop,kbot] are model interfaces (beginning at zero), 
! whereas in vd_lu_decomp they are expected as midpoints

! Compute rho at interfaces p/RT,  Ti = (Tm_k + Tm_k-1)/2,  interface k-
    do k = 2, pver
       do i = 1, ncol
          rho(i,k) = pi(i,k) * 2._r8 / (r * (t(i,k) + t(i,k-1)))
       end do
    end do
    do i = 1, ncol
       rho(i,pverp) = pi (i,pverp) / (r * t(i,pver))
    end do
    do i = 1, ncol
       rho(i,1) = pi (i,1) / (r * t(i,1))
    end do

! Calculate dt * (g*rho)^2/dp at interior interfaces,  interface k-
    tmpi2(:,:) = 0.0_r8
    do k = ktop+2 , kbot+1
       do i = 1, ncol
          tmpi2(i,k) = dt * (g * rho(i,k))**2 / (pm(i,k) - pm(i,k-1))
       end do
    end do

    do k = 1, pver
       do i = 1, ncol
          temp(i,k) =  dse(i,k)
       enddo
    enddo

! Decompose and solve the diffusion matrix
    call vd_lu_decomp(                                                &
         lchnk      ,ncol       ,                                     &
         zero       ,egwdffi    ,tmpi2       ,rdpm       ,dt         ,zero ,&
         ca         ,cc         ,term       ,ze         ,ktop+1      ,&
         kbot+1 )
    call vd_lu_solve(                                                 &
         lchnk      ,ncol       ,                                     &
         temp       ,ca         ,ze         ,term       ,             &
         ktop+1      ,kbot+1    ,zero   )

! Evaluated first part of tendency : diffusive term
   do k = 1, pver
      do i =1, ncol
         dttdf(i,k) = (temp(i,k) - dse(i,k)) / dt
      enddo
  enddo

! Evaluate second term: Conversion of kinetic energy into thermal
  do k = ktop+1, kbot
     do i = 1, ncol
        do l = -ngwv, ngwv
           dttke(i,k) = dttke(i,k) +  c(i,l) * gwut(i,k,l) 
        enddo
     enddo
  enddo

  do k=1, pver
     do i=1, ncol

        ttgw(i,k) = dttke(i,k) + dttdf(i,k)

     end do
  end do

  call outfld ('TTGWSDF', dttdf  / cpair, pcols, lchnk)
  call outfld ('TTGWSKE', dttke / cpair, pcols, lchnk)

  return
end subroutine gw_temp_tend


!==============================================================================


!+Mods for Beres04
subroutine gw_bgnd_beres (lchnk, ncol, &
     u, v, t, pm, pi, dpm, rdpm, piln,     &
     kldv, kldvmn, ksrc, ksrcmn, rdpldv, tau, ubi, ubm, xv, yv, &
     ngwv, zm, src_level)
!-----------------------------------------------------------------------
! Driver for multiple gravity wave drag parameterization.
!
! The parameterization is assumed to operate only where water vapor
! concentrations are negligible in determining the density.
!
! Beres, J.H., M.J. Alexander, and J.R. Holton, 2004: "A method of
! specifying the gravity wave spectrum above convection based on latent
! heating properties and background wind". J. Atmos. Sci., Vol 61, No. 3,
! pp. 324-337.
!
!-----------------------------------------------------------------------
  use phys_grid,    only : get_rlat_all_p, get_rlon_all_p
  use phys_grid,    only : get_lat_all_p, get_lon_all_p
  use phys_buffer,  only : pbuf

!------------------------------Arguments--------------------------------
  integer, intent(in) :: lchnk                  ! chunk identifier
  integer, intent(in) :: ncol                   ! number of atmospheric columns

  integer,  intent(in) :: ngwv                  ! number of gravity waves to use
  real(r8), intent(in) :: u(pcols,pver)         ! midpoint zonal wind
  real(r8), intent(in) :: v(pcols,pver)         ! midpoint meridional wind
  real(r8), intent(in) :: t(pcols,pver)         ! midpoint temperatures
  real(r8), intent(in) :: pm(pcols,pver)        ! midpoint pressures
  real(r8), intent(in) :: pi(pcols,0:pver)      ! interface pressures
  real(r8), intent(in) :: dpm(pcols,pver)       ! midpoint delta p (pi(k)-pi(k-1))
  real(r8), intent(in) :: rdpm(pcols,pver)      ! 1. / (pi(k)-pi(k-1))
  real(r8), intent(in) :: piln(pcols,0:pver)    ! ln(interface pressures)
  integer, intent(out) :: kldv(pcols)           ! top interface of low level stress divergence region
  integer, intent(out) :: kldvmn                ! min value of kldv
  integer, intent(out) :: ksrc(pcols)           ! index of top interface of source region
  integer, intent(out) :: ksrcmn                ! min value of ksrc

  real(r8), intent(in) :: rdpldv(pcols)        ! 1/dp across low level divergence region
  real(r8), intent(out) :: tau(pcols,-pgwv:pgwv,0:pver)! wave Reynolds stress
  real(r8), intent(out) :: ubi(pcols,0:pver)    ! projection of wind at interfaces
  real(r8), intent(out) :: ubm(pcols,pver)      ! projection of wind at midpoints
  real(r8), intent(out) :: xv(pcols)            ! unit vectors of source wind (x)
  real(r8), intent(out) :: yv(pcols)            ! unit vectors of source wind (y)
  real(r8), intent(in) ::  zm(pcols,pver)        ! midpoint heights
! C.-C. Chen
  integer, intent(out) :: src_level(pcols)      ! index of source (top of heating)
!---------------------------Local storage-------------------------------
  integer :: i,k,l                              ! loop indexes

  real(r8) :: rlon(pcols)                       ! longitude in radians for columns
  real(r8) :: rlat(pcols)                       ! latitude in radians for columns
  real(r8) :: tauback(pcols)                    ! background stress at c=0
  real(r8) :: u700(pcols)                       ! u at 700mb
  real(r8) :: v700(pcols)                       ! v at 700mb
  real(r8) :: al0                               ! Used in lat dependence of GW spec.
  real(r8) :: dlat0                             ! Used in lat dependence of GW spec.
  real(r8) :: flat_gw                           ! The actual lat dependence of GW spec.
  real(r8) :: pi_g                              ! 3.14........

!-------------------------- More local storage required for the Beres scheme -----------
  real(r8) :: netdt(pcols,pver)                 ! heating rate
  integer  :: lons(pcols)
  integer  :: lats(pcols)
  real(r8) :: curq(pver)           ! heating rate in given column
  real(r8) :: q0                                ! maximum heating rate
  real(r8) :: curu(pver)                        ! mean wind in a given column
                                                  ! projected in direction of 700 mb wind
  real(r8) :: depth                             ! heating depth
  real(r8) :: hdepth(pcols)                     ! heating depth array
  real(r8) :: maxq0(pcols)                      ! maximum heating rate array
!  real(r8) :: BSEMF(pcols)                      ! source total Eastward MF
!  real(r8) :: BSWMF(pcols)                      ! source total Westward MF

  real(r8) :: htop                              ! heating top
  integer  :: mini,maxi                         ! min and max heating range index
  integer  :: Umini,Umaxi                       ! Critical level filteringrange
  integer  :: maxsi                               ! maximum heating depth height index
  real(r8) :: uh                                ! mean wind in heating region
  real(r8) :: Umin, Umax                        ! Min and Max values of U in column
  real(r8) :: tau0(-PGWV:PGWV)                  ! source level tau for one column
  integer  CS                                   ! speed of convective cells relative to storm
  integer  shift                                ! how much to shift spectra to ground relative


  real(r8) :: ustorm                             ! storm speed
  real(r8), parameter :: CF = 20._r8                ! Heating rate conversion factor

  real(r8), parameter :: AL = 1.0e5_r8               ! Averaging length

  real(r8) :: dummy(pcols)
  real(r8) :: ccn
  character(len=80) fnum,fmt ,dumc1          ! dummies
  integer :: m
  real(r8), pointer, dimension(:,:) :: zmdt

  call get_lon_all_p(lchnk, ncol, lons)
  call get_lat_all_p(lchnk, ncol, lats)

  call get_rlat_all_p(lchnk, ncol, rlat)
  call get_rlon_all_p(lchnk, ncol, rlon)

!       do i = 1, ncol

!          write(iulog,*) 'ncol: ',ncol
!          write(iulog,*) 'Current, lat & lon: ',lats(i), lons (i)
!          write(iulog,*) 'Current, lat & lon: ',rlat(i)*180./3.14, rlon (i)*180./3.14

!       end do

  pi_g=4._r8*atan(1._r8)
  al0=40._r8*pi_g/180._r8
  dlat0=10._r8*pi_g/180._r8

  !----------------------------------------------------------------------
  ! Initialize tau array
  !----------------------------------------------------------------------

!  BSEMF(:) = 0.0
!  BSWMF(:) = 0.0
  tau(:,:,:)=0.0_r8
  netdt(:,:) = 0.0_r8
  hdepth(:) = 0.0_r8
  maxq0(:) = 0.0_r8
  dummy(:) = 0.0_r8
  tau0(:)  = 0.0_r8

  !---------------------------------------------------------------------------
  ! Determine 700 mb layer wind and unit vectors, then project winds.
  !---------------------------------------------------------------------------

  ! Just use the 700 mb interface values for the source
  ! wind speed and direction (unit vector).

  k = k700
  do i = 1, ncol
     u700(i) = u(i,k)
     v700(i) = v(i,k)
     ubi(i,k700) = sqrt (u700(i)**2 + v700(i)**2)
     xv(i) = u700(i) / ubi(i,k)
     yv(i) = v700(i) / ubi(i,k)
  end do

  ! Project the local wind at midpoints onto the direction of 700mb wind.
  do k = 1, pver
     do i = 1, ncol
        ubm(i,k) = u(i,k) * xv(i) + v(i,k) * yv(i)
     end do
  end do


  ! Compute the interface wind projection by averaging the midpoint winds.
  ! Use the top level wind at the top interface.
  ! Use the bottom level wind at surface

  do i = 1, ncol
     ubi(i,0) = ubm(i,1)
  end do
  do k = 1, pver-1
     do i = 1, ncol
        ubi(i,k) = 0.5_r8 * (ubm(i,k) + ubm(i,k+1))
     end do
  end do
  do i = 1, ncol
     ubi(i,pver) = ubm(i,pver)
  end do

  !---------------------------------------------------------------------------
  ! Find the index of z = 20 km (maximum possible heating range)
  ! using just the first column value of zm
  !---------------------------------------------------------------------------
  k=pver
  do
     if (zm(1,k).ge. 20000._r8) EXIT
     k=k-1
  end do
  maxsi = k


  !---------------------------------------------------------------------------
  ! Set up heating
  !---------------------------------------------------------------------------
  zmdt  => pbuf(idx_zmdt) %fld_ptr(1,:,:,lchnk,1)
  netdt(:pcols,:pver) = zmdt(:pcols,:pver)

  call outfld ('NETDT', netdt, pcols, lchnk)


  !-----------------------------------------------------------------------
  ! Gravity wave sources
  !!-----------------------------------------------------------------------
  !! Start LOOP over ALL COLUMNS
  !!-----------------------------------------------------------------------
  do i=1,ncol

     curq(1:pver) = netdt(i,1:pver)
     curu(1:pver) = ubm(i,1:pver)

     !-----------------------------------------------------------------------
     ! calculate HEATING DEPTH
     !
     ! heating depth is defined as the first  height range from the bottom
     ! in which heating rate is continuously positive
     !-----------------------------------------------------------------------
     k=pver


     ! Find first level of positive heating rate
     do
        if (k.le. maxsi) EXIT
        if (curq(k) .gt. 0.0_r8) EXIT
        k=k-1
     end do

     mini=k

     ! Find first level at which heating rate stops being positive
     do
        if (k.le. maxsi) EXIT
        if (curq(k) .le. 0.0_r8) EXIT
        k=k-1
     end do

     maxi=k
     depth = (zm(i,maxi)-zm(i,mini))/1000._r8
     htop  = zm(i,maxi)
     ! Depth is now in km

     hdepth(i) = depth
     ! Maximum heating rate
     q0 = MAXVAL(curq(maxi:mini))

     !output max heating rate in K/day
     maxq0(i) = q0*24._r8*3600._r8

     !-----------------------------------------------------------------------
     !Look up spectrum only if depth ge 2.5 km, else set tau0 = 0
     !-----------------------------------------------------------------------

     if ((depth .ge. 2.5_r8).AND.(ABS(rlat(i)).lt.(pi_g/2._r8))) then

        !-----------------------------------------------------------------------
        ! calculate MEAN WIND IN HEATING REGION relative to STORM CELLS
        !-----------------------------------------------------------------------

        ustorm = ubm(i,k700)               ! storm speed


        ! Calculate cell speed (CS) for ustrom > 10 m/s

        if (ABS(ustorm) .gt. 10._r8) then
           CS = SIGN(ABS(ustorm)-10._r8,ustorm)
        else
           CS = 0
        end if

        uh = 0.0_r8

        do k=maxi,mini
           uh = uh + ubm(i,k)/(mini-maxi+1)
        end do
        uh = uh - CS

        !-----------------------------------------------------------------------
        ! LOOK UP the spectrum using depth and uh
        ! Round both to the NEAREST WHOLE NUMBER
        ! SPECTRUM IS RElATIVE TO HEATING CELLS
        !-----------------------------------------------------------------------


        if(NINT(depth).gt.MAXH) then
           depth=MAXH*1.0_r8
           write(iulog,*) 'DEPTH exceeding table value'
        end if

        if(abs(NINT(uh)).gt.MAXUH) then
           uh=MAXUH*1.0_r8
           write(iulog,*) 'UH exceeding table value'
        end if

            tau0(:)=mfcc(NINT(depth),NINT(uh),:)

        ! NOW SHIFT SPECTRUM SO IT's RELATIVE TO THE GROUND

        shift = -NINT(CS/dc)

        tau0=CSHIFT(tau0,shift)

        ! Adjust Magnitude

        q0 = q0 * CF                                  ! multiply by conversion factor

        tau0(:)=tau0(:)*q0*q0/AL

        ! Adjust Magnitude for latitudinal dependence


!!$        flat_gw= 0.5*(1.+tanh((rlat(i)-al0)/dlat0)) + 0.5*(1.+tanh(-(rlat(i)+al0)/dlat0))
!!$        tau0(:)=(1-flat_gw)*tau0(:)


        ! Adjust for critical level filtering

!        Umin = Minval(ubm(i,bkbotbg:mini))            ! bkbotbg should be maxi
!        Umax = Maxval(ubm(i,bkbotbg:mini))            ! if launching waves from
                                                      ! top of heating

        Umin = Minval(ubm(i,maxi:mini))         
        Umax = Maxval(ubm(i,maxi:mini))

        Umini = max(NINT(Umin/dc),-PGWV)
        Umaxi = min(NINT(Umax/dc),PGWV)

        if (Umaxi.gt.Umini) then
           do l=Umini,Umaxi
              tau0(l) = 0.0_r8
           end do
        end if


     else

        tau0(:)=0.0_r8

     end if                                      !END IF depth ge 2.5


     ! Set GW source to the top of the heating
     !   write(iulog,*) 'maxi: ',maxi

!     ksrc(i) = bkbotbg
!     kldv(i) = bkbotbg
! C.-C. Chen
     ksrc(i) = maxi
     kldv(i) = maxi
     src_level(i) = maxi

     !launch waves at the same level for now

!     tau(i,-ngwv:ngwv,bkbotbg) = tau0(-ngwv:ngwv)
     tau(i,-ngwv:ngwv,maxi) = tau0(-ngwv:ngwv)

! Diagnostic variables
!     do l=1,ngwv
!        if (xv(i) .gt.0) then
!           BSEMF(i) = BSEMF(i)+ABS(tau0(l)*xv(i))
!           BSWMF(i) = BSWMF(i)+ABS(tau0(-l)*xv(i))
!        else
!           BSWMF(i) = BSWMF(i)+ABS(tau0(l)*xv(i))
!           BSEMF(i) = BSEMF(i)+ABS(tau0(-l)*xv(i))
!        end if

!     end do



  enddo

  !!-----------------------------------------------------------------------
  !! END LOOP over ALL COLUMNS
  !!-----------------------------------------------------------------------

  call outfld ('HDEPTH', hdepth, pcols, lchnk)
  call outfld ('MAXQ0', maxq0, pcols, lchnk)
!    call outfld ('BSEMF', BSEMF, pcols, lchnk)
!    call outfld ('BSWMF', BSWMF, pcols, lchnk)


  ! Determine the min value of kldv and ksrc for limiting later loops
  ! and the pressure at the top interface of the low level stress divergence
  ! region.

  ksrcmn = pver
  kldvmn = pver

  return
end  subroutine gw_bgnd_beres
!-Mods for Beres04


!==============================================================================
  subroutine vd_lu_decomp(                                          &
       lchnk      ,ncol       ,                                     &
       ksrf       ,kv         ,tmpi       ,rpdel      ,ztodt      , cc_top, &
       ca         ,cc         ,dnom       ,ze         ,ntop       , &
       nbot       )
!-----------------------------------------------------------------------
! Determine superdiagonal (ca(k)) and subdiagonal (cc(k)) coeffs of the 
! tridiagonal diffusion matrix. 
! The diagonal elements (1+ca(k)+cc(k)) are not required by the solver.
! Also determine ze factor and denominator for ze and zf (see solver).
!------------------------------Arguments--------------------------------
    integer, intent(in) :: lchnk                   ! chunk identifier
    integer, intent(in) :: ncol                    ! number of atmospheric columns
    integer, intent(in) :: ntop                    ! top level to operate on
    integer, intent(in) :: nbot                    ! bottom level to operate on

    real(r8), intent(in) :: ksrf(pcols)            ! surface "drag" coefficient
    real(r8), intent(in) :: kv(pcols,pverp)        ! vertical diffusion coefficients
    real(r8), intent(in) :: tmpi(pcols,pverp)      ! dt*(g/R)**2/dp*pi(k+1)/(.5*(tm(k+1)+tm(k))**2
    real(r8), intent(in) :: rpdel(pcols,pver)      ! 1./pdel  (thickness bet interfaces)
    real(r8), intent(in) :: ztodt                  ! 2 delta-t
    real(r8), intent(in) :: cc_top(pcols)          ! lower diagonal on top interface (for fixed ubc only)

    real(r8), intent(out) :: ca(pcols,pver)        ! -upper diagonal
    real(r8), intent(out) :: cc(pcols,pver)        ! -lower diagonal
    real(r8), intent(out) :: dnom(pcols,pver)      ! 1./(1. + ca(k) + cc(k) - cc(k)*ze(k-1))
                                                   ! 1./(b(k) - c(k)*e(k-1))
    real(r8), intent(out) :: ze(pcols,pver)        ! term in tri-diag. matrix system

!---------------------------Local workspace-----------------------------
    integer :: i                                   ! longitude index
    integer :: k                                   ! vertical index

!-----------------------------------------------------------------------
! Determine superdiagonal (ca(k)) and subdiagonal (cc(k)) coeffs of the 
! tridiagonal diffusion matrix. The diagonal elements  (cb=1+ca+cc) are
! a combination of ca and cc; they are not required by the solver.

    do k = nbot-1, ntop, -1
       do i = 1, ncol
          ca(i,k  ) = kv(i,k+1)*tmpi(i,k+1)*rpdel(i,k  )
          cc(i,k+1) = kv(i,k+1)*tmpi(i,k+1)*rpdel(i,k+1)
       end do
    end do

! The bottom element of the upper diagonal (ca) is zero (element not used).
! The subdiagonal (cc) is not needed in the solver.

    do i=1,ncol
       ca(i,nbot) = 0._r8
    end do

! Calculate e(k).  This term is 
! required in solution of tridiagonal matrix defined by implicit diffusion eqn.

    do i = 1, ncol
       dnom(i,nbot) = 1._r8/(1._r8 + cc(i,nbot) + ksrf(i)*ztodt*gravit*rpdel(i,nbot))
       ze(i,nbot)   = cc(i,nbot)*dnom(i,nbot)
    end do
    do k = nbot-1, ntop+1, -1
       do i=1,ncol
          dnom(i,k) = 1._r8/ (1._r8 + ca(i,k) + cc(i,k) - ca(i,k)*ze(i,k+1))
          ze(i,k) = cc(i,k)*dnom(i,k)
       end do
    end do
    do i=1,ncol
       dnom(i,ntop) = 1._r8/ (1._r8 + ca(i,ntop) + cc_top(i) - ca(i,ntop)*ze(i,ntop+1))
    end do

    return
  end subroutine vd_lu_decomp

!==============================================================================
  subroutine vd_lu_solve(                                           &
       lchnk      ,ncol       ,                                     &
       q          ,ca         ,ze         ,dnom       ,             &
       ntop       ,nbot       ,cd_top )
!-----------------------------------------------------------------------
! Solve the implicit vertical diffusion equation with zero flux boundary conditions.
! Actual surface fluxes are explicit (rather than implicit) and are applied separately. 
! Procedure for solution of the implicit equation follows Richtmyer and 
! Morton (1967,pp 198-200).

! The equation solved is
!
! -ca(k)*q(k+1) + cb(k)*q(k) - cc(k)*q(k-1) = d(k), 
!
! where d(k) is the input profile and q(k) is the output profile
!
! The solution has the form
!
! q(k)  = ze(k)*q(k-1) + zf(k)
!
! ze(k) = cc(k) * dnom(k)
!
! zf(k) = [d(k) + ca(k)*zf(k+1)] * dnom(k)
!
! dnom(k) = 1/[cb(k) - ca(k)*ze(k+1)] =  1/[1 + ca(k) + cc(k) - ca(k)*ze(k+1)]

! Note that the same routine is used for temperature, momentum and tracers,
! and that input variables are replaced.
!------------------------------Arguments--------------------------------
    integer, intent(in) :: lchnk                   ! chunk identifier
    integer, intent(in) :: ncol                    ! number of atmospheric columns
    integer, intent(in) :: ntop                    ! top level to operate on
    integer, intent(in) :: nbot                    ! bottom level to operate on

    real(r8), intent(in) :: ca(pcols,pver)         ! -upper diag coeff.of tri-diag matrix
    real(r8), intent(in) :: ze(pcols,pver)         ! term in tri-diag solution
    real(r8), intent(in) :: dnom(pcols,pver)       ! 1./(1. + ca(k) + cc(k) - ca(k)*ze(k+1))
    real(r8), intent(in) :: cd_top(pcols)          ! cc_top * ubc value

    real(r8), intent(inout) :: q(pcols,pver)       ! constituent field

!---------------------------Local workspace-----------------------------
    real(r8) :: zf(pcols,pver)                     ! term in tri-diag solution

    integer i,k                                    ! longitude, vertical indices

!-----------------------------------------------------------------------

! Calculate zf(k).  Terms zf(k) and ze(k) are required in solution of 
! tridiagonal matrix defined by implicit diffusion eqn.
! Note that only levels ntop through nbot need be solved for.

    do i = 1, ncol
       zf(i,nbot) = q(i,nbot)*dnom(i,nbot)
    end do
    do k = nbot-1, ntop+1, -1
       do i=1,ncol
          zf(i,k) = (q(i,k) + ca(i,k)*zf(i,k+1))*dnom(i,k)
       end do
    end do

! Include boundary condition on top element
    k = ntop
    do i=1, ncol
       zf(i,k) = (q(i,k) + cd_top(i) + ca(i,k)*zf(i,k+1))*dnom(i,k)
    end do

! Perform back substitution

    do i=1,ncol
       q(i,ntop) = zf(i,ntop)
    end do
    do k=ntop+1, nbot, +1
       do i=1,ncol
          q(i,k) = zf(i,k) + ze(i,k)*q(i,k-1)
       end do
    end do

    return
  end subroutine vd_lu_solve

!===============================================================================
end module gw_drag
